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Laminar  and  Turbulent  Flows  Past  Two 
Dimensional  and  Axisymmetric  Bodies 

ABSTRACT 

A  numerical  study  of  laminar  and  turbulent  flows  past  two  dimensional  bodies  and 
axisymmetric  bodies  is  presented.  Numerical  methods  are  developed  to  solve  Navier- 
Stokes  equations  for  two  dimensional  and  axisymmetric  flows  in  the  arbitrary  geometries. 
The  complex  physical  geometry  is  resolved  by  use  of  numerically  generated,  body-fitted 
coordinates.  The  governing  equations  are  written  in  the  transformed  domain  using  the 
orthogonal  velocity  components  as  dependent  variables  for  momentum  equations.  The 
governing  equations  are  discretized  using  both  the  finite  analytic  method  and  the  finite 
volume  method.  Both  one  velocity  staggered  grid  method  and  two  velocities  staggered 
method  are  employed  for  grid  arrangements.  The  velocity  and  pressure  coupling 
techniques  in  these  grid  arrangements  are  presented.  The  solution  procedure  of  the 
SIMPLER  numerical  algorithm  is  used  with  a  parabolic  marching  technique  and  a  global 
pressure  calculation  method.  For  turbulent  flow  calculations,  both  the  k-|  turbulence 
model  and  the  two-layer  turbulence  model  are  used.  (wo 

Calculations  are  performed  for  laminar  and  turbulent  flows  past  a  finite  flat  plate  and 
turbulent  flow  past  axisymmetric  bodies  with  different  solution  domains,  numerical 
methods  and  turbulence  models.  Calculations  include  the  development  of  a  wake  function 
method  for  the  prediction  of  turbulent  wake  of  a  flat  plate,  prediction  %  of  ^  aminar  and 
turbulent  flows  past  a  finite  flat  plate,  predictions  of  tu  'bulent  flow  past  axisymmetric 
bodies  by  the  wall  function  method  and  by  the  two-layer  turbulence  model  and  predictions 
of  turbulent  flow  past  finite  axisymmetric  bodies.  Comparisons  of  predictions  by  finite 
analytic  method  with  those  by  finite  volume  method  are  made  for  some  calculations.  All 
the  predicted  results  were  compared  with  available  experiment  measurements.  Good 
agreements  between  the  predicted  results  and  measured  data  were  obtained. 
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PREFACE 


This  monograph  entitled  "The  Finite  Analytic  Method  and  Its  Applications"  contains  the 
continuing  developments  of  the  numerical  method  called  the  "finite  analytic"  method.  The  finite 
analytic  method  was  developed  in  early  1977  when  I  and  Dr.  Peter  Li,  Aen  a  graduate  student,  had 
difficulty  in  obtaining  a  converged  solution  from  a  system  of  finite  difference  ^gelnaic  equations 
derived  for  the  Navier-Stokes  equations  for  two-dimensional  turbulent  flow  with  a  second  order 
turbulence  noodeL  I  conceived  the  finite  analytic  method  one  night  and  solved  the  simple  two- 
dimensional  Laplace  equation.  Dr.  Li  then  carried  the  finite  analytic  method  to  the  unsteady 
diffusion  equation  and  nonlinear  ordinary  differential  equation.  The  basic  idea  of  the  finite  analytic 
method  is  to  derive  the  ^proximate  algebraic  representation  of  a  governing  linear  or  nonlinear 
differential  equation  firom  the  local  analytic  solution.  The  local  analytic  solution  is  obtained  for  a 
small  element  of  the  total  solution  domain  in  which  the  governing  equation,  if  nonlinear,  is 
linearized.  The  local  analytic  solution  is  then  expressed  in  an  algebraic  form.  The  system  of  local 
analytic  algebraic  equations  is  solved  to  provide  the  numerical  solution  of  the  problem. 

Subsequently,  in  a  span  of  ten  years,  many  students  took  interest  in  the  development  of  the 
finite  analytic  method  and  its  ^plications  to  fluid  mechanics  and  heat  transfer  problems.  They  are 
Hamid  Naseri-Neshat,  Hamn-Ching  Chen,  Kuo-San  Ho,  Young  Hwan  Yoon,  Che  Hsi  Yu,  Sen- 
Ming  Chang,  Wenchung  Chen,  Zahed  Mohammad  Sheil^leslami,  Tzong-Shyan  Wung,  Wu  Sun 
Chen,  Kemakolam  Obasih,  Seok  Ki  Choi,  Vahid  Talaie,  Ramiro  Humberto  Bravo,  and  Hakan 
Aksoy.  They  all  contribute  to  the  further  development  of  the  finite  analytic  method.  In  1983  a 
monogr^h  entitled  "The  Finite  Analytic  Method"  was  published  as  part  of  the  Iowa  Institute  of 
Hydraulic  Research  Report.  There  are  a  total  of  seven  volumes  in  the  monograph  entitled  "The 
Hnite  Analytic  Method."  The  special  feature  of  the  finite  analytic  method  is  that  when  it  is  applied 
to  solve  the  Navier-Stokes  equations,  it  provides  relatively  accurate  representations  of  the 
convection  term.  The  fiinite  analytic  method  ntn  tmly  provides  accurate  simulation  of  the 
upwinding  effea  because  of  its  analytic  nature  of  the  solutions,  but  it  also  provides  automatic, 
gradual  shifits  of  the  upwinding.  Con^uently,  the  finite  anal}^c  numeric^  solutions  of 
convection  diffusion  equations  minimize  the  false  numerical  diffusion  that  would  occur  in  the 
upwinding  difference  used  in  other  numerical  methods  and  provide  stable  and  fast  convergency  of 
the  solutions. 

Six  years  have  passed  since  the  publication  of  the  monograph  entitled  'Tinite  Analytic 
Method."  During  these  years  the  emphasis  has  been  to  apply  the  finite  analytic  method  to  various 
problems  and  to  develop  a  more  fiiendly  code  fw  users.  This  monograph  Aus  is  developed  to 
document  the  efforts  m^e  in  the  applications  of  the  finite  analytic  method  along  with  its  further 
improvements. 

I  would  like  to  thank  my  colleagues,  V.C.  Patel,  K.B.  Chandran,  T.F.  Smith,  A.  T.  Chwang, 
F.  Stem,  and  K.  Atkinson  for  their  interest  in  the  development  and  application  of  the  finite  analytic 
method.  I  would  like  to  acknowledge  the  support  of  Drs.  William  D.  McNally,  Peter  M  Sockol, 
Gary  Johnson,  and  JJ.  Adamezyk  of  NASA  Lewis  Center  for  their  support  in  the  early 
development  of  the  finite  analytic  method.  My  thanks  also  go  to  Dr.  Oscar  P.  Manley  of  the  U.S. 
Department  of  Energy  for  his  supi^  of  the  application  of  the  finite  analytic  method  in  energy 
related  problems.  Recent  applications  of  the  finite  analytic  method  to  ship  hydrodynamics 
problems  have  been  support^  by  the  Navy  GHR  Grant  No.  N0014-84-0068,  by  the  Naval  Sea 
S3rstem  Command  GHR  Grant  No.  N(X)168-86-J-(X)19  administered  by  the  David  Taylor  Naval 
Ship  Research  ^  Development  Center,  and  by  a  recent  grant  from  the  Office  of  Navd  Research 
under  the  Applied  Hydrodynamics  Research  Grant  No.  N00167-86-K-(X)19. 


In  Vdume  I  (tf  the  Hnite  Analytic  Method  and  Its  Application,  laminar  and  turbulent  flows 
past  a  two-dimensional  and  axisymmetiic  body  are  studied  I  would  like  to  thank  Dr.  Seok  Ki 
Qkh  for  applying  the  finite  analytic  method  to  con:q)lex  flow  problems.  Without  his  participatitm, 
the  finite  ai^ytic  method  would  not  have  been  developed  and  understood  as  it  is  today.  I  sbdl  be 
luqtpy  to  receive  any  discussion  or  criticism  on  the  finite  analytic  metiiod  I  attempt  to  prewnt  the 
finite  analytic  metii^  and  its  applications  in  a  systematic  documentation.  Whatever  possible,  all 
the  listings  of  the  codes  (fevelc^)^  associated  witii  the  finite  analytic  method  are  oiclosed  so  tiiat 
readers  may  apply  the  finite  analytic  method  in  their  problems. 
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ABSTRACT 


A  numerical  study  of  laminar  and  turbulent  flows  past 
two  dimensional  bodies  and  axisymmetric  bodies  is  presented. 
Numerical  methods  are  developed  to  solve  Navier-Stokes 
equations  for  two  dimensional  and  axisymmetric  flows  in  the 
arbitrary  geometries.  The  complex  physical  geometry  is 
resolved  by  use  of  numerically  generated,  body-fitted 
coordinates .  The  governing  equations  are  written  in  the 
transformed  domain  using  the  orthogonal  velocity  components 
as  dependent  variables  for  momentum  equations.  The  governing 
equations  are  discretized  using  both  the  finite  analytic 
method  and  the  finite  volume  method.  Both  one  velocity 
staggered  grid  method  and  two  velocities  staggered  method  are 
employed  for  grid  arrangements .  The  velocity  and  pressure 
coupling  techniques  in  these  grid  arrangements  are  presented. 
The  solution  procedure  of  the  SIMPLER  numerical  algorithm  is 
used  with  a  parabolic  marching  technique  and  a  global 
pressure  calculation  method.  For  turbulent  flow 
calculations,  both  the  k  -  e  turbulence  model  and  the  two- 
layer  turbulence  model  are  used. 

Calculations  are  performed  for  laminar  and  turbulent 
flows  past  a  finite  flat  plate  and  turbulent  flow  past 

iii 


axisymmetric  bodies  with  different  solution  domains, 
numerical  methods  and  turbulence  models.  Calculations 
include  the  development  of  a  wake  function  method  for  the 
prediction  of  turbulent  wake  of  a  flat  plate,  predictions  of 
laminar  and  turbulent  flows  past  a  finite  flat  plate, 
predictions  of  turbulent  flow  past  axisymmetric  bodies  by  the 
wall  function  method  and  by  the  two-layer  turbulence  model 
and  predictions  of  turbulent  flow  past  finite  axisymmetric 
bodies.  Comparisions  of  predictions  by  finite  analytic 
method  with  those  by  finite  volume  method  are  made  for  some 
calculations.  All  the  predicted  results  were  compared  with 
available  experimental  measurements.  Good  agreements  between 
the  predicted  results  and  measured  data  were  obtained. 
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CHAPTER  I 
INTRODUCTION 

I . 1  Purpose  of  Study 

Recent  demands  in  advanced  ship  design  have  prompted 
many  researchers  to  study  flow  past  ship  propulsion 
systems .  Since  most  ship  propellers  are  located  at  the 
downstream  end  of  ship  body,  the  study  of  flow  past  ship 
bodies  is  important  for  the  design  of  optimal  ship  shape 
and  propellers.  However,  the  present  state-of-the-art  of 
prediction  methods  of  flow  past  complex  geometries  is  not 
yet  capable  of  caculating  the  entire  flow  field  past  ship 
bodies.  The  investigations  of  flows  past  axisymmetric 
bodies  can  provide  a  fundamental  understanding  of  flow 
field  past  more  complex  ship  bodies.  Studies  of  fluid 
flows  past  two  dimensional  bodies  such  as  hydrofoils  and 
appendages  are  also  important  for  ship  design  since  the 
understanding  of  flow  field  past  appendages  will  provide 
more  realistic  simulation  of  inlet  flow  conditions  for  the 
propeller. 

The  purpose  of  the  present  study  is  to  develop 
numerical  calculation  methods  for  predicting  laminar  and 
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turbulent  flows  past  two  dimensional  and  axisymmetric 
bodies .  The  present  study  starts  with  the  prediction  of 
laminar  and  turbulent  flows  past  a  thin  flat  plate  and 
then  turbulent  flows  past  axisymmetric  bodies  of  DTNSRDC 
(David  Taylor  Naval  Ship  Research  and  Development  Center) 
are  investigated  in  detail.  These  geometries  are  chosen 
not  only  because  their  geometries  approximately  resemble 
the  appendages  or  ship  body  but  also  because  of  the 
availability  of  extensive  experimental  data  to  compare 
with  the  predicted  results.  The  shapes  of  these  bodies 
are  presented  with  generated  numerical  grid  in  the 
following  chapter. 

In  the  following  sections,  the  theoretical  backgroud 
of  present  investigation  is  presented.  It  includes  the 
review  of  numerical  methods,  grid  generation  techniques, 
calculation  methods  for  incompressible  flow  in  arbitrary 
geometry,  and  turbulence  model. 

_ Theoretical  Background 

Fluid  flows  past  two  dimensional  or  axisymmetric 
bodies  are  also  encountered  in  a  variety  of  practical 
engineering  applications,  for  example,  aerodynamics  and 
turbomachineries.  The  most  reliable  information  may  be 
obtained  by  direct  measurements.  However,  experimental 


3 


investigations  are  in  most  case  expensive  and  time 
consuming.  There  also  exist  difficulties  of  measurement 
in  many  situations,  and  measuring  instruments  are  not  free 
from  errors . 

On  the  other  hand,  with  the  rapid  development  of 
computer  technology  in  computational  speed  and  size  of 
storage,  the  computational  methods  became  practical  to 
simulate  various  flow  phenomena  encountered  in  many 
engineering  problems.  During  the  past  twenty  years,  the 
computational  methods  for  incompressible  flow  have  been 
rapidly  developed  accompanied  with  the  advances  of 
practical  turbulence  model.  Various  computational  methods 
have  been  proposed,  tested  and  refined  to  a  stage  which 
may  significantly  impact  the  design  of  many  engineering 
problems .  Engineers  and  researchers  begin  to  realize  that 
the  CFD  (Computational  Fluid  Dynamics)  is  perhaps  a  cost 
effective  and  convenient  way  of  analyzing  the  complex 
engineering  problems . 

Two  numerical  methodologies  are  the  most  commonly 
employed  and  developed,  namely  finite  difference  method 
and  finite  element  method.  Most  flows  in  the  engineering 
problems  have  rather  arbitrary  geometries .  When  the 
practical  usage  of  calculation  methods  in  arbitrary 
geometry  is  considered,  the  finite  element  method  appears 


4 


to  be  the  natural  choice.  However,  the  finite  element 
method  still  has  unresolved  difficulties  in  the  simulation 
of  fluid  flow  phenomena  although  this  method  is  widely 
used  in  other  areas  such  as  solid  mechanics.  By  its 
rather  simplicity  and  efficiency,  the  finite  difference 
method  has  been  the  most  widely  used  in  the  computational 
fluid  dynamics.  Among  the  various  finite  difference 
methods,  the  method  developed  by  Patankar  and  Spalding 
[1,2]  seems  to  be  the  most  popular  in  the  calculation  of 
incompressible  flows. 

In  addition  to  the  finite  difference  method  and  the 
finite  element  method,  a  finite  analytic  method  has  been 
developed  by  Chen  and  Chen  [3] .  This  method  has  been 
tested  for  the  several  practical  engineering  problems 
[4,5,6]  and  proved  to  be  accurate  and  stable.  In  the 
finite  analytic  method,  the  algebraic  representation  of  a 
certain  nodal  value  with  neighbouring  nodal  values  are 
obtained  by  a  local  analytic  solution  of  the  governing 
equation.  The  novel  nature  of  this  method  is  that  the 
numerical  false  diffusion  problem  in  solving  the  Navier- 
Stokes  equations,  which  is  considered  as  the  most  serious 
problem  for  nearly  all  the  numerical  schemes  of  finite 
difference  method,  is  satisfactorily  minimized  by  using 
all  the  surrounding  nodal  points.  The  well  developed 
calculation  procedure  for  the  finite  difference  method  in 
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solving  the  Naviers-Stokes  equations  can  be  readily 
applied  to  finite  analytic  method  without  any 
difficulties.  Both  the  finite  analytic  method  and  finite 
difference  method  are  used  for  the  several  calculations 
presented  in  this  study.  Comparisions  of  predictions  by 
these  two  methods  are  given  for  certain  problems . 

Recently,  several  calculation  methods  have  been 
developed  that  employ  boundary-fitted,  curvilinear 
coordinate  to  appropriately  handle  fluid  problems  in 
complex  geometries .  The  technique  of  generating  an 
appropriate  boundary  fitted  curvilinear  coordinate  and  the 
development  of  related  calculation  procedure  associated 
with  using  the  curvilinear  coordinate  system  have  been 
subjects  of  numerous  studies.  Although  the  calculation 
methods  for  incompressible  flow  in  the  orthogonal  grid 
system  have  been  well  developed  and  refined,  those  for 
non-orthogonal  grid  system  seem  to  be  not  yet  fully 
established.  Detailed  discussion  of  currently  available 
calculation  methods  for  incompressible  flow  associated 
with  using  the  curvilinear  coordinte  system  will  be  given 
in  the  following  section. 

Due  to  its  importance  in  the  engineering  application, 
several  numerical  and  experimental  studies  of  laminar  and 
turbulent  flows  past  two  dimensional  or  axisymmetric 
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bodies  have  been  made  during  the  last  decade.  Since  this 
problem  includes  the  leading  edge  interaction,  boundary 
layer  development  on  the  body,  the  trailing  edge 
interaction  and  wake  development,  it  has  been  subject  of 
testing  of  several  calculation  methods  and  turbulence 
models.  At  present,  several  reliable  experimental 
informations  are  available  which  can  be  used  to  verify  the 
accuracy  of  the  calculation  schemes  or  turbulence  models. 
Detailed  discussion  of  previous  numerical  and  experimental 
studies  is  given  in  the  chapter  where  each  specific 
problem  is  mentioned. 

Previous  numerical  and  experimental  studies  indicate 
that  although  the  first  order  boundary  layer  theory 
adequately  describes  flow  fields  over  the  middle  of  body, 
it  fails  to  describe  flow  fields  near  the  leading  and 
trailing  edge  of  a  body  and  near  the  tail  region  of 
axisymmetric  bodies.  Two  calculation  methods  have  been 
widely  used  in  order  to  appropriately  handle  these  flow 
situations.  The  first  method  is  the  viscous-inviscid 
interaction  method  in  which  boundary  layer  solutions  are 
matched  with  the  external  inviscid  flows  by  interactive 
means.  The  second  approach  is  the  global  numerical 
solutions  of  elliptic  or  partially  parabolic  form  of 
Navier-Stokes  equations.  The  latter  approach  is  more 
attractive  since  it  does  not  require  another  numerical 
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solution  of  inviscid  flow  and  is  adopted  in  the  present 
study. 

In  many  previous  numerical  studies  of  flow  past  ship 
bodies,  calculations  are  started  either  at  the  middle  of 
the  body  or  at  the  trailing  edge  of  the  body  with  upstream 
conditions  specified  by  well  documented  profiles  from 
boundary  layer  theory  or  experimental  data.  Since  these 
calculations  do  not  include  the  leading  edge  interaction 
and  initial  development  of  boundary  layer  over  the  body, 
the  resulting  solutions  are  dependent  on  the  upstream 
conditions  or  the  location  of  upstream  solution  boundary. 

A  more  reliable  numerical  solution  may  be  obtained  by 
including  the  leading  edge  within  the  solution  domain  with 
prescribing  uniform  flow  conditions  at  the  inlet  and  with 
the  solution  domain  large  enough  to  capture  the  whole 
viscous-inviscid  interaction.  However,  inclusion  of 
leading  edge  in  the  solution  domain  may  cause  the 
difficulties  of  turbulence  modelling  in  the  transition 
region  for  turbulent  flow  calculation.  In  the  present 
study,  the  upstream  solution  boundary  is  made  to  locate 
either  at  the  uniform  flow  region  far  upstream  of  body  or 
at  the  middle  of  the  body  according  to  the  objective  of 
the  problem  considered. 
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Since  the  problems  considered  in  the  present  study 
have  rather  arbitrary  geometries  which  require  the  use  of 
boundary-fitted  curvilinear  coordinate  system,  a  brief 
review  of  methods  for  the  generation  of  boundary-fitted 
coordinates  and  available  calculation  methods  for 
incompressible  flow  in  arbitrary  geometryis  is  presented. 

A  brief  discussion  of  available  turbulence  models  will  be 
also  presented. 

1.3  Generation  of  Boundary-Fitted  Coordinahea 

A  boundary-fitted  coordinate  system  is  defined  as  a 
curvilinear  coordinate  system  in  which  the  boundaries  of 
the  physical  domain  concide  with  a  curvilinear  coordinate 
line  or  surface.  The  generation  of  an  appropriate 
computational  grid  is  very  important  in  order  to  achieve 
the  accuracy  and  better  convergence  of  numerical  solution. 
Boundary-fitted  curvilinear  coordinates  can  be  generated 
orthogonally  or  non-orthogonally .  The  orthogonal  grid  has 
the  advantages  that  the  transformed  partial  differential 
equations  are  simpler  as  the  non-orthogonal  terms  do  not 
appear.  This  advantage  results  in  fewer  computing  time, 
fast  convergence  as  well  as  better  stability  and  accuracy 
of  the  solution.  However,  orthogonal  coordinates  are 
often  difficult  to  generate  and  it  is  sometimes  impossible 
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to  generate  the  orthogonal  coordinate  for  the  three 
dimensional  situations.  An  additional  disadvantage  for 
the  generation  of  orthogonal  coordinate  is  the  fact  that 
there  is  no  control  on  the  location  of  the  grid  points 
along  the  boundary  which  may  deteriorate  the  grid  spacing 
near  th^  solid  boundary  where  the  appropriate  grid  spacing 
is  required.  Non-orthogonal  coordinates  are  much  easier 
to  generate  and  can  be  produced  also  for  three  dimensional 
situations.  Further,  concentration  of  grid  lines  in  the 
regions  where  good  resolution  is  required  is  much  easier 
to  achieve.  However,  non-orthogonal  coordinate  introduce 
additional  terms  such  as  cross  derivatives  which  may 
reduce  both  convergence  and  stability  of  numerical 
solution  and,  thereby,  increase  the  computing  time. 

There  exist  several  ways  of  the  generating  boundary- 
fitted  coordinates.  The  available  grid  generation 
techniques  can  be  broadly  classified  into  two  categories, 
algebraic  methods  and  partial  differential  equation 
methods.  The  multi-surface  technique  of  Eiseman  [7]  is  a 
typical  method  for  algebraic  method  in  which  simple 
functions  are  combined  to  generate  grids  for  complex 
geometries.  In  this  approach,  the  grid  is  generated  by 
joining  corresponding  points  on  the  inner  and  outer 
boundaries  by  polynomial  curves.  In  general,  the 
algebraic  grid  generation  techniques  are  very  attractive 
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by  the  reasons  that  they  do  not  require  a  numerical 
solution  of  partial  differential  equations  and  a  precise 
control  over  the  grid  spacing  is  possible.  However,  when 
this  method  is  applied  to  a  complex  geometry,  non¬ 
smoothing  grids  may  be  produced. 

Boundary-fitted  coordinates  can  also  be  generated  by 
the  solution  of  partial  differential  equations  in  which 
the  dependent  variables  are  the  grid  coordinates  in  the 
physical  domain.  Conformal  mapping  is  a  special  case  of 
grid  generation  technique  using  the  Laplace  equation. 

Much  of  the  works  have  been  done  by  Thompson  and  his 
coworkers  [8] .  The  spacing  between  the  grid  lines  is 
strongly  dependent  on  the  equations  being  solved  and 
weakly  dependent  on  the  boundary  point  distributions .  The 
equations  can  be  easily  changed  by  varying  the  values  of 
the  source  term  or  the  grid  control  functions  in  the 
Poisson  equations.  These  grid  control  functions  can  be 
selected  to  appropriately  attract  the  coordinate  line 
toward  particular  grid  points  or  grid  lines  or  repel  the 
coordinate  lines.  Thompson  et  al.  [8]  have  suggested  a 
series  of  exponential  functions  as  source  terms.  However, 
the  selection  of  grid  control  function  is  c[uite  dependent 
on  the  problem  being  solved. 
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In  the  present  study,  all  the  numerical  grids  are 
generated  by  the  solution  of  the  Poisson  type  of  partial 
differential  equations  with  grid  control  functions 
specified  properly  depending  on  the  physical  situations 
being  solved.  Details  of  grid  generation  techniques 
adopted  in  the  present  study  are  given  in  chapter-II. 

1.4  Calculation  Methods  for  Incompressible  Flow  in 

Complex  Geometries 

1.4.1  Calculation  Methods  for  Incompressible  Flow  in 
Orthogonal  Coordinate  System 

The  computational  method  for  the  prediction  of  fluid 
flow  in  the  orthogonal  grid  system  is  well  developed  and 
have  been  widely  used  for  a  variety  of  problems .  The 
orthogonal  grid  based  calculation  method  is  not  much 
depart  from  that  for  Cartesian  counterparts.  The  well 
developed  calculation  procedure  for  Cartesian  coordinate 
system  can  be  easily  modified  with  a  proper  specification 
of  geometric  coefficients  and  source  terms.  Most  of 
previous  calculations  use  the  staggered  grid  arrangement 
and  the  velocity  components  normal  to  the  control  volume 
surfaces  (contravariant  velocity  components)  are  selected 
as  the  dependent  variables  in  the  momentum  equations. 
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Pope  [9]  used  such  a  orthogonal  based  calculation 
procedure  to  compute  turbulent  recirculating  flows  in 
diffusers,  Gosman  and  Rapley  [10]  applied  a  similiar 
method  to  the  calculation  of  fully  developed  flow  in  the 
ducts  of  arbitrary  cross  section.  Recently,  Raithby  et 
al.  [11]  have  presented  a  new  calculation  procedure  in 
which  the  stresses  are  retained  in  the  discretization 
equation  instead  of  the  usual  procedure  of  substituting 
them  in  terms  of  strain  rates . 

However,  the  orthogonal  coordinate  based  methods  have 
limitations  in  the  application  due  to  the  difficulty  of 
generating  appropriate  numerical  grid.  For  the  sake  of 
general  application,  the  calculation  methods  for  non- 
orthogonal  coordinate  system  should  be  investigated. 

1.4.2  Calculation  Methods  for  Incompressible  Flow  in 
Non-orthogonal  Coordinate  System 

Unlike  the  orthogonal  coordinate  based  methods, 
several  calculation  methods  for  the  non-orthogonal 
coordinate  system  have  been  proposed  during  the  late  1970s 
and  during  the  1980s.  The  differences  among  these  methods 
are  the  choice  of  dependent  variables  in  the  momentum 
equations,  the  grid  arrangements  and  the  treatment  of 
pressure  and  velocity  coupling.  One  interesting  fact  is 
that  none  of  these  methods  have  superiority  over  other 


methods  to  justify  its  general  usage.  Each  of  the  methods 
has  its  own  merits  and  demerits  which  should  be  judged  in 
terms  of  the  problem  being  considered.  As  compared  with 
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the  method  based  on  the  orthogonal  coordinate  system, 
these  methods  are  still  in  a  developing  stage.  Among  the 
several  calculation  methods  based  on  non-orthogoal 
coordinate  system,  five  typical  methods  which  use 
different  grid  arrangements  or  dependent  varibles  as  shown 
in  Fig.I-1  are  discussed  below. 

Non-Staagered  Grid  Method 

In  the  first  grid  arrangement,  which  is  shown  in  the 
Fig.I-l-a,  all  the  dependent  variables  are  computed  and 
stored  at  the  same  location  and  the  Cartesian  velocity 
components  are  used  as  dependent  variables  in  the  momentum 
equation.  This  method  has  been  named  differently 
according  to  the  author's  choice.  For  example,  P  -  Q 
scheme  by  Hsu  [12];  Momentum  Interpolation  Method  by 
Majumdar  [13];  Pressure  Weighted  Interpolation  Method  by 
Miller  and  Schmidt  [14];  Collocation  Method  by  Peric  et 
al.  [15];  and,  Non-Staggered  Grid  Method  by  Rhie  and  Chow 
[16] .  Since  this  method  uses  Cartesian  velocity 
components  as  dependent  variables  for  the  momentum 
equation,  the  governing  equation  is  simple  and  the 
curvature  terms  which  are  extremely  grid  dependent  can  be 


avoided.  Moreover,  a  strongly  conservative  form  of 
transport  equation  is  possible  which  is  a  desired 
characteristics  for  finite  volume  method.  As  is  discussed 
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in  Patankar  [2],  this  grid  arrangement  may  give  rise  to  a 
checkerboard  pressure  oscillation.  To  overcome  this 
checkerboard  splitting  of  pressure  field,  a  special  kind 
of  interpolation  method  for  evaluating  velocity  components 
at  the  control  volume  cell  surface  has  been  devised.  All 
the  transport  equations  are  solved  at  the  center  point  of 
control  volume  cell  and  the  cell  surface  velocities  for 
the  derivation  of  pressure  or  pressure  correction  equation 
are  obtained  by  a  linear  interpolation  using  neighbouring 
nodal  values  except  the  pressure  gradient  term  in  the 
momentum  equation  is  evaluated  by  the  method  employed  in 
the  staggered  grid  approach  to  prevent  the  splitting  of 
the  pressure  field.  Thus,  this  method  indirectly  uses  the 
staggering  idea  and  it  was  questioned  that  this  method  can 
be  named  as  non-staggered  method  [17] .  Since  all  the 
variables  are  calculated  at  the  same  point,  the 
coefficients  of  algebraic  equations  are  the  same  for  all 
the  variables  which  may  reduce  the  computer  storage  and 
computing  time,  especially  for  the  three  dimensional 
situations.  This  feature  will  lead  to  more  cost  effective 
calculation  when  the  multigrid  method  is  employed.  The 
earlier  method  by  Rhie  and  Chow  [16]  or  Peric  [18]  has 
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been  blamed  for  the  fact  that  the  resulting  converged 
solution  is  quite  relaxation  factor  or  time  step 
dependent.  Recently,  Majumdar  [13,19]  and  Miller  and 
Schmidt  [14]  removed  this  problem  by  a  rigorous 
derivation.  This  method  has  been  successfuly  applied  to  a 
variety  of  engineering  problems  by  many  researchers. 
However,  Miller  and  Schmidt  [14]  find  that  this  method  may 
result  in  a  physically  unrealistic  solution  in  regions 
where  there  is  a  strong  pressure  gradient .  This  method 
also  suffers  from  the  lack  of  diagonal  dominance  of 
pressure  correction  equation  when  it  is  applied  to 
strongly  non-orthogonal  grid.  This  method  has  been 
applied  to  a  variety  of  problems;  Rhie  and  Chow  [16],  Rhie 
and  coworkers  [20,21],  Peric  [18],  Han  [22],  Majumdar 
[13,19],  Majumdar  et  al.  [17],  Rodi  et  al .  [23],  Miller 
and  Schmidt  [14]  and  Peric  et  al.  [15]  .  An  excellent 
discussion  of  this  method  is  given  in  Majumdar  et  al . 

[17]  . 

One  Velocity  Staggered  Grid  Method 

Shyy  et  al.  [24]  developed  a  calculation  method 
employing  the  conventional  staggered  grid  usually  adopted 
in  the  Cartesian  coordinate  system.  As  shown  in  the 
Fig.I-l-b,  the  dependent  variables  in  the  momentum 
equations  are  the  Cartesian  velocity  components  and  only 
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one  velocity  component  is  stored  at  the  each  cell  surface. 
Since  only  one  velocity  component  is  available  at  each 
cell  surface,  the  other  unavailable  velocity  component 
which  is  necessary  to  compute  the  mass  flow  rate  is 
usually  obtained  through  the  linear  interpolation  of 
neighbouring  velocity  components  which  are  obtained  in  the 
previous  iteration.  In  addition  to  this  defficiency, 
further  approximations  have  to  be  made,  for  example 
neglect ion  of  non-orthogonal  pressure  gradient  terms,  to 
provide  the  diagonal  dominant  pressure  correction  equation 
in  strongly  non-othorgonal  grid.  A  more  serious  drawback 
of  this  method  is  that  if  the  flow  direction  is  normal  to 
the  direction  of  velocity  component  stored  at  the  cell 
surface,  zero  convective  mass  fluxes  across  the  control 
volume  faces  may  occur  and  the  pressure  gradient  driving 
the  velocity  components  will  be  completely  wrong,  which 
may  cause  undesired  pressure  splitting.  However,  use  of 
this  method  is  justified  when  the  flow  direction  is  not 
much  deviated  from  the  direction  of  velocity  components 
stored  at  the  cell  surface.  Due  to  the  simplicity  of  this 
method,  many  researchers  used  this  method  for  the  problem 
to  which  this  method  can  be  applied.  For  example,  Braaten 
and  Shyy  [251,  Nakayama  [26],  Ha  [27],  Chen  and  Patel  [4], 
Cheng  and  Chen  [5],  Obashi  and  Chen  [6],  Patel  et  al .  [28] 
and  Choi  and  Chen  [29] . 


Discussions  of  the  previous  method  suggest  that  in 


order  to  calculate  flow  fields  around  arbitrary 
geometries,  both  of  the  Cartesian  velocity  components 
should  be  stored  at  each  cell  surface  as  shown  in  Fig.I-1- 
c,  while  pressure  and  turbulent  quantities  are  stored  at 
the  center  point  of  control  volume  cell.  Maliska  and 
Raithby  [30]  adopted  this  configuration  in  the  development 
of  three  dimensional  parabolic  calculation  procedure  for 
the  fluid  flow  in  the  arbitrary  cross  sectioned  ducts. 

This  method  is  also  not  free  from  difficulties.  In  order 
to  obtain  the  diagonal  dominant  pressure  correction 
equation  in  the  strongly  non-orthogonal  grid,  the  non- 
orthogonal  pressure  gradient  term  should  be  neglected. 
Storing  the  two  velocity  components  at  each  control  volume 
surface  requires  a  complex  programing  and  more  computer 
storage  and  more  computing  time.  These  disadvantages 
became  more  grave  in  three  dimensional  situations.  This 
method  is  not  commonly  used  due  to  forementioned 
disadvantages . 

Algebraic  Manipulation  Method 

An  alternative  to  the  Cartesian  velocity  components 
is  to  use  the  curvilinear  velocity  components  as  the 
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dependent  variables  in  the  momentum  equations.  The 
contravariant  or  the  covariant  velocity  components  will  be 
the  natural  choices  due  to  their  relationship  with  the 
grid  lines.  Karki  and  Patankar  [31]  developed  a 
calculation  method  using  covariant  velocity  components  as 
the  dependent  variables  in  the  momentum  equations .  The 
discretization  equations  and  source  terms  for  covariant 
velocity  components  were  obtained  by  the  algebraic 
manipulation  of  discretization  equation  and  source  terms 
based  on  the  Cartesian  velocity  components.  As  shown  in 
the  Fig.I-l-e,  each  covariant  velocity  component  is  stored 
at  the  control  volume  cell  surface  using  staggering  idea. 
The  usage  of  covariant  velocity  components  as  the 
dependent  variable  in  the  momentum  equation  enables  the 
strongly  diagonal  dominant  pressure  correction  equation 
even  in  the  strongly  non-orthogonal  grid  situations.  The 
complicated  tensor  algebra  is  avoided  through  a  novel 
algebraic  manipulation  of  discretization  equation  and  the 
resulting  algebraic  equations  seem  to  be  cost  effective 
and  less  grid  dependent .  The  success  of  this  method  is 
largely  dependent  on  the  accuracy  and  grid  dependency  of 
the  algebraic  manipulation  of  source  term  of  momentum 
equation  which  is  considered  as  the  most  serious  problem 
when  the  contravariant  or  covariant  velocity  components 
are  adopted  as  a  dependent  variables.  This  method  should 
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be  tested  for  complicated  three  dimensional  situations  to 
justify  its  general  usage.  This  method  has  been 
successfully  applied  to  calculation  of  compressible  and 
incompressible  in  two  dimensional  complex  geometries 
[32,33].  Detailed  study  of  this  method  is  provided  by 
Karki  [34] . 

It  is  worthy  to  mention  the  calculation  method  based 
on  the  full  transformation  of  dependent  variables  for 
momentum  equations  using  a  complicated  tensor  algebra. 
Demirdzic  et  al.  [35]  developed  such  a  calculation  method 
using  contravariant  velocity  components  as  dependent 
variables  for  momentum  equations.  Fig.I-l-d  shows  that 
the  only  one  contravariant  velocity  component  is  stored  at 
the  each  control  volume  cell  surface  and  the  staggering 
idea  is  adopted.  This  method  requires  a  complicated 
tensor  algebra  and  more  storages  for  the  geometric 
coefficients.  Moreover,  the  resulting  governing  equation 
is  extremely  complicate  and  the  numerical  solution  is 
quite  grid  dependent.  This  fact  restricts  the  usage  of 
this  method  in  certain  situations.  Like  all  the  methods 
discussed  previously,  this  method  also  gives  a  nine  point 
pressure  correction  equation  which  lacks  diagonal 
dominance  in  strongly  non-orthogonal  grid  situations  if 
the  non-orthogonal  pressure  gradient  terms  are  not 
neglected.  The  advantage  of  smaller  numerical  diffusion 
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in  this  approach  must  be  judged  against  the  serious 
disadvantages  mentioned  above.  This  method  has  been  used 
by  few  people.  For  example,  Demirdzic  et  al.  [35], 
Demirdzic  [36],  Gal-Chen  et  al.  [37],  Stern  et  al.  [38] 
and  Richmond  [39] .  An  excellent  study  of  this  method  is 
given  by  Demirdzic  [36] . 

Methods  of  Vanka  et  al.  [40]  and  Faghri  et  al.  [41] 
deserve  to  receive  attention.  However,  these  methods  are 
also  not  free  from  many  difficulties  mentioned  in  the 
previous  methods  and  are  not  commonly  used  by  many 
researchers . 

The  one-velocity  staggered  grid  method  is  adopted  for 
the  most  of  calculations  of  the  present  study  since  this 
method  is  very  easy  to  program  and  the  flow  direction  is 
not  much  deviated  from  the  direction  of  velocity  component 
stored  at  cell  surface  with  an  appropriate  grid 
generation.  The  two-velocities  staggered  grid  method  is 
also  used  for  some  calculations  in  the  present  study  with 
some  modifications  made  from  the  original  method  of 
Maliska  and  Raithby  [30] . 

I . 5  Turbulence  Models 

During  the  last  decade,  many  turbulence  models 
ranging  from  the  Prandtl  mixing-length  model  [42]  to  the  k 
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-  e  turbulence  model  [43]  were  used  for  a  variety  of 
practical  calculations  and  tested  in  order  to  find  out 
their  range  of  applicability.  Some  progress  has  been 
made,  but  to  date  no  model  has  been  found  to  be  both 
accurate  and  general.  The  models  widely  used  in 
applications  have  been  based  on  the  Boussinesq  (1877) 
assumption  which  relates  the  turbulent  stresses  to  mean 
rates  of  strain  by 

au.  au.s  2 

<  aq  a^  )  -  3  ^ij 

where  is  turbulent  eddy  viscosity  and  k  is 
turbulent  kinetic  energy. 

For  a  suitable  charateristic  length  scale  1^  and 
velocity  scale  v^,  dimensional  analysis  suggests  that 
turbulent  eddy  viscosity  can  be  evaluated  as 

(1-2) 

The  difference  among  the  existing  models  based  on  the 
Boussinesq  assumption  is  how  velocity  scale  v^  and  length 
scale  1^  are  specified.  Experimental  evidence  indicates 
that  the  Boussinesq  assumption  is  reasonably  valid  in  many 
flow  circumtances  with  exception  of  some  flow  situations. 
Models  based  on  the  Boussinesq  assumption  are  frequently 
called  eddy  viscosity  model. 
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Turbulence  models  are  often  classified  according  to 
the  number  of  supplementary  partial  differential  equations 
that  must  be  solved  in  order  to  supply  the  modelling 
parameters .  Since  details  of  these  models  are  available 
in  the  various  literatures,  for  example  Chen  [44],  only  a 
brief  discussion  of  current  available  turbulence  models 
will  be  given. 

1.5.1  Zero  Equation  Turbulence  Model 

In  the  zero  equation  model  or  algebraic  model,  such 
as  proposed  by  Cebeci  and  Smith  [45]  or  Baldwin  and  Lomax 
[46],  the  eddy  viscosity  is  prescribed  from  length  scale 
and  velocity  gradient. 

The  length  scale  1^  is  usually  obtained  from  the 
empirical  correlations  developed  for  the  boundary  layers. 
The  algebraic  models  have  proven  to  be  accurate  and 
reliable  for  relatively  simple  flows.  The  disadvantage  of 
this  model  is  the  lack  of  generality  that  the  length  scale 
which  is  needed  for  specification  of  eddy  viscosity  is 
function  of  various  parameters  and  it  is  not  easy  to 
specify  an  appropriate  length  scale  over  the  whole  flow 
field,  especially  for  the  separated  flow  region  or  the 
wake  region.  However,  these  models  can  be  employed  all 
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the  way  to  the  wall  where  the  viscosity  effect  is 
important  and  the  numerical  solutions  of  partial 
differential  equations  are  not  needed  for  the 
specification  of  length  scale. 

1.5.2  One  Equation  Turbulence  Model 

This  model  is  based  on  the  suggestion  of  Prandtl  and 
Kolomogorov  that  the  turbulent  velocity  scale  is 
proportional  to  the  square  root  of  turbulent  kinetic 
energy  k,  implying; 

Vt  =  k^^^  1^  (1-4) 

This  model  requires  the  numerical  solution  of 
turbulent  kinetic  energy  equation  which  can  be  derived 
from  the  Navier-Stokes  equations.  The  turbulent  kinetic 
energy  equation  contains  a  diffusion  term,  a  generation 
term  and  a  dissipation  terra  which  must  be  modelled  through 
additional  assumptions  such  as; 


which  requires  specification  of  another  length  scale 

There  exist  several  one  equation  models  [47,48,49] 
and  the  differences  among  the  these  models  are  mainly  how 
the  length  scales  1^,  Ig  are  specified.  In  general,  the 
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predictions  of  the  one  eqpjation  model  have  been  marginally 
better  than  those  obtained  from  zero  equation  models . 


1.5.3  k  -  e  Turbulence  Model 

The  k  -  e  turbulence  model  [43]  is  usually  called  two 
equation  model.  In  this  model,  the  eddy  viscosity  is 
related  to  the  turbulent  kinetic  energy  k  and  its 
dissipation  rate  e  ; 

2 


''t  “ 

^  M-  e 


(1-6) 


The  turbulent  dissipation  rate  e  is  assumed  to  be 
related  to  the  length  scale  through 
.3/2 


^t  “ 


(1-7) 


The  turbulent  kinetic  energy  k  and  its  dissipation 
rate  e  are  determined  from  two  modelled  equations .  Thus 
the  eddy  viscosity  and  the  length  scale  are  obtained  from 
the  numerical  solutions  of  partial  differential  equations. 
Since  this  model  does  not  require  algebraic  specification 
of  length  scale,  it  can  be  applied  to  several  flow 
situations . 


The  k  -  e  turbulence  model  can  only  be  applied  to  the 
fully  turbulent  region  and  calculations  of  the  viscosity- 
affected  near  wall  region  is  avoided  by  the  use  of  wall 
functions.  The  wall  function  method  [43]  has  been  found 
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to  be  quite  successful  for  attached  flow,  it  is  not 
suitable  for  some  flow  situations.  The  wall  functions 
relating  the  velocity  and  turbulent  quantities  at  the 
second  grid  point  to  friction  velocity  heavily  lean  on  the 
assumptions  of  a  logarithmic  velocity  distribution  and  the 
validity  of  local  equilibrium  of  turbulence.  These 
assumptions  are  not  valid  in  the  separated  flow  region. 
Furthermore,  the  extension  of  wall  function  method  to 
unsteady  flow  and  three  dimensional  flow  requires 
additional  assumptions  which  are  not  substantiated  by 
experiment  or  theory  [50] .  An  alternative  to  the  use  of 
wall  function  method  will  be  to  employ  the  turbulence 
model  which  are  valid  all  the  way  to  the  wall. 

1.5.4  Low  Reynolds  Number  Turbulence  Model 

Several  low  Reynolds  number  extensions  of  the  k  -  6 
turbulence  model  [51,52,53]  which  can  be  applicable  to  the 
viscous  affected  layer  have  been  proposed.  These  models 
introduce  damping  functions  in  order  to  achieve  the 
observed  reduction  of  turbulent  transport  very  near  the 
wall.  However,  these  models  have  an  undesirable  feature 
of  requiring  very  high  numerical  resolution  near  the  wall 
region  because  of  the  steep  gradient  of  the  dissipation 
rate  e.  It  was  found  that  all  these  models  perform  rather 
poorly  in  adverse  pressure  gradient  boundary  layers  [54] 
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and  in  separated  flows  [55].  Patel  et  al.  [50]  reviewed 
the  several  available  low  Reynolds  number  k  -  e  turbulence 
models  and  concluded  that  none  of  the  these  models  can  be 
used  confidently  even  in  the  boundary  layer  calculations. 

1.5.5  Two-Layer  Turbulence  Model 

In  order  to  reduce  the  computational  effort  and  to 
preserve  numerical  stability,  the  well  established  length 
scale  distribution  near  the  wall  region  can  be  introduced 
to  the  calculation  of  the  vicosity-affected  near  wall 
region  instead  of  using  the  wall  function  method  while  the 
outer,  fully  turbulent  region  is  resolved  by  use  of  the  k 

-  e  turbulence  model.  This  model  is  named  as  two-layer 
model  and  it  is  usually  used  with  the  combination  of  the  k 

-  e  turbulence  model  and  zero  equation  or  one  equation 
turbulence  model.  Some  successes  in  the  prediction  of 
complex  turbulent  flows  by  use  of  this  model  is  reported 
[56,57,58,59,60]  . 

In  the  present  study,  the  k  -  £  turbulence  model  with 
wall  function  method  will  be  employed  for  most  of 
turbulent  flow  calculations  while  the  two-layer  model  is 
applied  to  some  flow  situations  involving  the  pressure 
gradient  and  body  curvature  effects  to  improve  the 
accuracy  of  the  calculations. 
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I • 6 _ Outline  of  Study 

In  chapter  II,  mathematical  and  numerical 
formulations  used  in  the  present  study  are  presented.  The 
governing  equations,  boundary  conditions,  methods  of 
generation  of  boundary-fitted  coordinates  and 
dicretization  schemes  are  discussed. 

In  chapter  III,  numerical  aspects  are  discussed. 

These  include  the  coupling  between  the  continuity  equation 
and  momentum  equations,  numerical  algorithm  and  overall 
solution  procedures. 

Chapters  IV  to  VII  are  devoted  to  the  numerical 
calculations  which  include  development  of  the  wake 
function  method  for  the  prediction  of  flat  plate  wake, 
prediction  of  laminar  and  turbulent  flows  past  a  finite 
flat  plate,  numerical  solutions  of  turbulent  flow  past 
axisymmetric  bodies  by  wall  function  method  and  by  two 
layer  model,  prediction  of  turbulent  flow  past  finite 
axisymmetric  bodies.  Some  comparisions  of  the  predictions 
by  finite  analytic  method  with  those  by  finite  volume 
method  are  presented.  All  the  calculted  results  are 
compared  with  available  experimental  data  to  demonstrate 
the  applicability  and  accuracy  of  the  numerical  scheme. 
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In  chapter  VIII,  a  brief  sunimary  of  the  present  study 
is  given  and  some  suggestions  for  future  work  are 
discussed. 


CHAPTER  II 


FORMULATION  OF  PROBLEM 

II  ..1 _ Introduction 

In  this  chapter,  the  mathematical  and  numerical 
formulation  of  the  present  study  is  presented.  This 
chapter  starts  with  the  mathematical  formulation  of  the 
problem  which  includes  governing  equations,  boundary 
conditions  and  solution  domains.  Then  the  grid  generation 
techniques  and  discretization  methods  employed  in  the 
present  study  are  presented. 

The  solution  domains  are  made  large  enough  to  capture 
the  whole  important  feature  of  flow  field.  The  numerical 
grids  are  generated  through  the  numerical  solutions  of 
Poisson  equations  in  which  the  dependent  variables  are  the 
physical  coordinate  of  solution  domain.  The  governing 
equations  are  written  in  the  transformed  domain  using  the 
partial  transformation  where  the  original  Cartesian  or 
cylindrical  velocity  components  are  left  as  dependent 
variables  for  momentum  equations.  These  transformed 
governing  equations  are  discretized  using  both  finite 
analytic  method  and  finite  volume  method. 


The  continuity  equation  and  general  form  of  transport 


equations  for  the  incompressible  flow  in  Cartesian  or  in 
cylindrical  coordinate  system  with  the  k  -  e  turbulence 
model  can  be  written  as  follow,' 

^(pru)+^(prv)=0  (II-l) 

^(pr(l)  )  +^(pru(j)  )  +^(prv(l)  ) 

a  ^  ^  a  T.  ^ 


In  these  equations,  ^  represents  the  general 
dependent  variables  (  u,  v,  k,  e  ) ,  p  is  the  density  of 
fluid,  is  the  effective  diffusion  coefficient,  (  u,  v 
are  the  velocity  components  in  (  x,  y  )  directions  and  S 
denotes  the  source  term  for  the  variable  (j) .  For  two 
dimensional  flow,  r  =  1  and  for  axisymmetric  flow,  y  =  r 
Details  of  the  definition  of  all  the  variables  are  given 


) 

4> 


in  Table  1 . 


I I. 3.1  Solution  Domains 


All  the  calculations  presented  in  the  present  study 
can  be  broadly  classified  into  two  cartegories;  half  body 
calculations  and  full  body  calculations.  The  difference 
between  these  two  calculations  is  the  location  of  the 
upstream  solution  boundary  as  shown  in  Fig.II-1.  In  the 
half  body  calculations,  calculations  start  at  the  middle 
of  body  with  inlet  conditions  specified  by  the  well 
documented  profiles  from  the  boundary  layer  theory.  In 
the  full  body  calculations,  calculations  start  far 
upstream  of  body  with  inlet  conditions  specified  by 
uniform  flow  conditions.  The  solution  domain  in  the 
normal  direction  and  in  the  downstream  direction  is  made 
large  enough  to  capture  the  whole  viscous-inviscid 
interaction  and  wake  development.  The  location  of 
upstream  inlet  solution  boundary  in  the  full  body 
calculations  is  also  made  far  away  from  the  body  t' 
capture  the  whole  leading  edge  interaction  The  details 
of  the  solution  domains  used  in  the  present  calculations 
are  given  in  Tables  2  to  6. 
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1 1. 3. 2  Boundary  Conditions 

The  boundary  conditions  for  all  the  problems 
considered  in  the  present  study  can  be  stated  as  follows; 

Inlet  or  Upstream  Boundary  {  x  =  ) 

The  upstream  inlet  solution  boundary  is  located 
either  at  the  middle  of  the  body  or  at  the  uniform  flow 
region  far  upstream  of  the  body.  The  u-velocity  component 
and  turbulent  quantities  k,  e  are  prescribed  from  the 
simple  flat-plate  correlations  for  the  half  body 
calculations  and  from  uniform  flow  conditions  for  the  full 
body  calculations.  For  v-velocity  component,  the 
condition  of  =  0  is  specified  for  all  the  calculations. 

Outlet  or  Downstream  Boundary  (  x  =  ) 

The  downstream  outlet  condition  is  placed  in  the  wake 
region  far  downstream  of  body.  The  zero  pJ^essure-gradient 
condition  p^  =  0  is  specified  at  the  outlet  boundary.  For 
the  exit  boundary  conditions  for  transport  quantities  (  u, 
V,  k,  e  )  ,  “  ^'^x  “  ^x  “  ^x  “  ®  conditions  are  specified. 

Body  Surface 


The  no-slip  boundary  condition  (  u  =  v  =  0  )  is 
specified  for  all  the  laminar  flow  calculations  and  for 
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the  turbulent  flow  calculations  when  the  two  layer  model 
is  employed  to  calculate  the  flow  field  all  the  way  to  the 
wall.  For  all  the  other  turbulent  flow  calculations,  the 
wall  function  method  is  employed  in  order  to  avoid  the  use 
of  large  number  of  grids  to  resolve  the  steep  gradient  in 
the  near  wall  region  and  uncertanties  of  the  turbulence 
model  in  the  viscous-affected  region. 

Upper  Boundary  (  y  =  y^  ) 

The  upper  boundary  is  placed  at  a  large  distance  from 
the  body  where  the  flow  is  assumed  to  be  uniform.  The 
uniform  boundary  condition  is  ; 

u  =  1,  p  =  0,  ky  =  Ey  =  0.  (II-3) 

The  v-velocity  component  in  this  boundary  is  obtained 
from  the  law  of  conservation  of  mass  within  the 
computational  control  volume  cell  during  the  solution 
process . 

Wake  Centerline  (  y  =  0  ) 

The  boundary  conditions  along  the  wake  centerline  are 
the  symmetry  conditions  and  are  specified  as 

=0,  Uy  =  ky  =  Ey  = 


V 


0. 


(II-4) 
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II.4.  Numerical  Grid  Generation 


The  boundary-fitted  coordinate  generation  technique 
proposed  by  Thompson  ec  al.  [8]  is  employed  for  the 
generation  of  numerical  grid  in  the  present  study.  The 
complex  solution  domain  in  physical  coordinates  (  x,  y  ) 
is  transformed  to  a  simple  rectangular  domain  in  non- 
orthogonal,  body-fitted  coordinates  (  Tj  )  by  the 
solutions  of  the  Poisson  equations; 


;)v2  ^  r  dv  ^ 


3x2  "  r  ay  3y2 

a^Tl  a^  ^  a2ij 

3x2  r  ay  3y2 


(II-5) 

(11-6) 


1  2 

where  f  and  f  are  the  grid  control  functions  which 
are  chosen  priori  to  obtain  the  desired  distribution  of 
grid  points  in  the  physical  (  x,  y  )  domain.  In  these 
equations,  r  =  1,  a^,  =  0  for  a  Cartesian  coordnate  system 
and  y  =  r  and  a^,  =  1  for  a  cylindrical  coordinate  system. 
The  origin  of  coordinate  system  is  located  at  the  leading 
edge  of  body. 


The  numerical  grid  is  generated  by  solution  of 
inverted  forms  of  these  equations  where  the  dependent 
variables  are  the  physical  coordinates  of  solution  domain. 


Il3£x  229£x 


2g 


i2_3£x_ 


1^  2ax 
94  9n 


0 


(II-7) 
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11^  0^2^  ^  12_a2j^  2^  ^ 

9^2  9t|2  9^9r|  9^  9ii  ^ 

where 


11 

2  2  2 
r  (  +  y^  ) 

(II-9) 

22 

2  2  2 
r  (  +  y^  ) 

(II-IO) 

12 

2 

-  r  (  XpX  +  y  y  ) 

(ii-ii) 

2 

■  r 

(  ) 

(11-12) 

The  grid  generation  techniques  outlined  below  are 
essentially  the  same  as  those  reported  in  Chen  and  Patel 
[4]  or  Cheng  and  Chen  [5].  However,  some  modifications 
are  made  for  certain  problems . 

II. 4.1  Generation  of  Numerical  Grids  for  Flat  Plate 

Non-uniform  rectangular  numerical  grids  are  generated 
for  the  calculation  of  laminar  and  turbulent  flow  past  a 
flat  plate  with  a  proper  grid  concentration  made  near  the 
wall  and  near  the  leading  and  trailing  edges.  The 
numerical  grid  in  the  normal  direction  is  generated 
algebraically  instead  of  solving  the  differential 
equation.  The  first  few  grids  near  the  plate  are 
generated  uniformly  and  next  the  other  grids  are  generated 
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non-uniformly  with  the  grid  expansion  ratio  about  1.2  to 
1.3.  The  size  of  the  first  grid  in  the  normal  direction 
from  the  plate  is  adjusted  according  to  the  physical 
situations  being  considered.  The  numerical  grid  in  x- 
direction  is  generated  by  the  numerical  solution  of  the 
following  differential  equation  which  is  obtained  from 
Eq. (II-7) . 


11  d^x  _1  dx  - 


(11-13) 


The  grid  control  function  f  is  selected  by  following 
equations  which  is  previously  used  by  Chen  and  Patel  [4] 
or  Cheng  and  Chen  [5]; 


f^/2g^^  =  Ai 


0.0  =  =  0.5 


sin(  TC  )  0.5  »  Z^  =  2 


(11-14) 


A2  sin(  1Z  Z2  )  0.0  =  Z2*1.5 


-  Ac 


1.5  =  Z« 


where 


Zi  = 


iiL 


Zo  = 


(11-15) 


The  terms  and  ^2  correspond  to  the  leading  and 
trailing  edge  stations  as  shown  in  Fig.II-2  and  A^^  and  A2 
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are  positive  constants  which  can  be  adjusted  to  yield  the 
desired  grid  concentration  near  the  leading  and  trailing 
edges . 

The  partial  view  of  generated  numerical  grids  used 
for  laminar  and  turbulent  flow  calculations  are  shown  in 
Fig.II-3.  As  shown  in  the  figures,  an  appropriate  grid 
concentration  was  made  near  the  wall  region  and  near  the 
leading  and  trailing  edge  regions  and  the  desired  grid 
expansion  in  the  wake  region  is  also  achieved. 


II. 4. 2  Generation  of  Numerical  Grids  for  Axisymmetric 

Bodies 


Body  fitted  coordinates  for  axisymmetric  bodies  are 
generated  by  the  numerical  solutions  of  Poisson  equations. 
Eg. (II~7)  and  Eq. (II-8) .  In  the  present  study,  the 
constant-x  lines  are  chosen  as  constant-^  lines  to 
facilitate  the  comparision  of  numerical  results  with 
experimental  data.  Then,  Eq. (I 1-7)  and  Eq. (I 1-8)  can  be 
rewritten  as  following  simplified  form. 


Il3£x  ,l9x 
Q  - TT  +  f  -  =  0 

lia2r  ^  2232r  ^  _  12a2r  idr  2dr 

«  W  ^  ^  rrr-  .  t  ..  .  F  . 

where  =  f^  - 


d^an  a^  an 


(11-16) 

(11-17) 


rr. 


(11-18) 
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The  x-coordinates  are  algebraically  generated  instead 
of  solving  Eq. (11-16) .  The  x-coordinate  over  the  body  and 
in  the  wake  region  near  the  trailing  edge  of  body  (  x/L  < 
1.1  )  is  generated  uniformly  while  x-coordinate  in  the 
wake  region  away  from  the  body  (  x/L  >1.1  )  is  generated 
nonuniformly  with  grid  expansion  ratio  about  1,2  to  1.3. 
The  grid  control  function  f^  is  then  determined  from 
following  equation  using  the  x-coordinate  already 
generated. 

.11^^ 


f  =  -  g 


(11-19) 


With  fl  specified  by  Eq. (11-19),  Eq. (11-17)  can  be 
solved  with  a  proper  specification  of  grid  control 
function  .  The  resulting  numerical  solution  gives  the 
grid  distribution  in  the  radial  direction. 


Along  the  inlet  and  outlet  boundary  where  the  grid  is 
orthogonal,  Eq. (11-17)  can  be  written  as  following. 


22d^r  odz  _ 


at  ^  “  i  ^  “  ^max  (11-20) 


Thus  the  grid  control  function  F  along  these 
boundaries  can  be  determined  from  Eq. (11-20)  if  the  r- 
coordinate  at  these  boundary  is  provided. 

r. 


Fa(t1) 


r2a,n) 


(11-21) 


(11-22) 


r 

FgCTl)  ”  ®  “  9  — 321 1  ^  ^  ^max 

r\ 

where  ^  -  1  is  the  upstream  inlet  station  and  ^  -  ^max 
is  the  downstream  outlet  station. 

2 

The  value  of  F  inside  the  solution  domain  can  be 
obtained  by  a  linear  combination  of  and  Fg. 


Fa(t1) 

A 

A 

> 

(11-23) 

Fc(^/'n) 

%K<%  <%B 

Fed!) 

^  ^max 

where 

Fc(^.Tl)  =  [(x(^b)-x(^))F^(T1)  +  (x(^)-x(^a))F3(T1)] 

/  tx(^B)-x(^A)  ]  (11-24) 

and  4a  determined  from  the  body  shape.  In 

the  case  of  half  body  calculation,  ^A  is  placed  at  the 
station  where  radius  of  body  is  initially  changed  and  4b 
is  placed  at  the  trailing  edge  of  body  as  shown  in  Fig. II- 
2.  However,  the  specification  of  F  (4^11)  is  quite  problem 
dependent  and  can  be  adjusted  to  obtain  the  desired  grid 
distribution  depending  on  the  shape  of  geometry. 

With  the  specification  of  proper  grid  control 
1  2 

functions  f  and  F  ,  Eq. (11-17)  is  solved  by  finite 
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difference  method.  Eq. (11-17)  is  discretized  as  following 
using  the  exponential  scheme  of  Spalding  [4] . 

11  22 

(  2a  g  coth  a  +  2b  g  coth  b  ) ^  ^  r^  ^ 

-  (  2a  each  a  >5,,  (  e'®  ) 

+  (  2b  9^^  csch  b  (  e'”  +  e*’  ) 

+  2  (9^2)^  ,,  (  ) 

(11-25) 


where 

2^--4:U.n 

g  g 

The  above  system  of  algebraic  equation  is  solved  by 
the  tridiagonal-matrix  algorithm. 

The  body  fitted  coordinate  system  for  full  body 

calculation  was  obtained  with  a  slightly  different 

specification  of  grid  control  function  F  (^,1))  .  In  order 

to  attract  grid  lines  toward  body  surface  in  the  leading 

edge  region  of  body,  the  following  grid  control  function 
2 

F  (4/'n)  was  devised. 

F^(4,ti)  -  [  (d(^B)-d(^)  (d(^)-d(^A)  )Fb(t1)  ] 


/  [d(^B)-d(^A)  1 


(11-27) 


41 


where  d(4)  -  -  r(^,l)  (11-28) 

In  Eq. (11-27),  the  station  is  the  station 

where  d(^)  is  minimum  and  the  station  ^  is  the  station 

where  d(^)  is  maximum  as  shown  in  Fig.II-2.  The  grid 

control  functions  at  these  stations,  ^^(71)^  Fg(T|)  are 

evaluated  using  Eq. (11-21)  and  Eq. (11-22)  as  the  same  way 

as  done  in  half  body  calculation.  The  normal  reference 

station  11  =  ^ref  shown  in  Fig.II-2  was  adjusted 

to  obtain  the  desired  grid  concentration  toward  body  and 

Tl  ^  =  T\  -  5  is  used  for  most  calculations. 

'ref  'max 

As  seen  from  the  Fig.II-4  and  Fig.II-5,  the  boundary- 
fitted  coordinates  generated  by  above  method  evolve 
smoothly  from  the  body  into  the  walce  and  the  desired  grid 
expansion  in  the  tail  region  and  in  the  wake  region  is 
obtained. 


II. 5.1  Partially  Parabolic  Form  of  Finite  Analytic  Method 
in  the  Cartesian  Coordinate  system 

The  non-dimensionized  form  of  general  transport 
equation  in  the  Cartesian  coordinate  system  can  be  written 


as  follows; 
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d<t)  d<|>  d<|>  1  92({)  32<j>  ^ 

~  +  u^  +  v^  -  ^  )  +  s'*’ 


dt 


dx 


ay 


3x2  9y2 


where  is  source  term  for  variable  (j)  and 
effective  Reynolds  number  for  variable  (|>  defined  as 


(11-29) 

is 


R  =  R  = 
u  V 


(1/Re  +  v^.) 


Ok 

(Ok/Re  +  v^.) 
(Oe/Re  +  v^) 


(11-30) 

(11-31) 

(11-32) 


If  we  locally  linearize  the  convective  velocity  in 
the  governing  equation  by  the  value  at  calculation  point  P 
as  shown  in  Fig.II-6  and  neglect  the  diffusion  of 
transport  quantities  in  the  flow  direction,  the  governing 
equation  can  be  approximately  written  in  a  following 
simplified  form: 

32<j>  3<|> 

^  -  2B  +  F  (11-33) 


where 


2B  -  (R^v)p 

-  <Vf  <  ♦p  -  1>w 


At 


(11-34) 


<Vp 

(11-35) 


and 


<V>P 


Ax 


(11-36) 
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-  (R4,)p  (11-37) 

In  these  equations,  the  subscript  P  denotes  the 
calculation  point  and  the  superscript  n-1  denotes  the 
previous  time  step  value .  The  unsteady  term  and  the  x- 
derivative  of  <j>  are  approximated  by  backward  and  upwind 
finite  difference  respectively.  Eq. (11-30)  is  solved 
analytically  with  boundary  conditions  specified  by 
following  equations. 

<j»(h)  =  (Hjj,  <{»(-k)  =  (l)g  (11-39) 

Evaluation  of  <})  at  calculation  point  P  gives 
following  algebraic  equations. 

fp  *  »N  ♦n  *  *3  - 


where 


^  exp(2Bh)  -  1 
exp(2Bh)  -  exp(-2Bk) 


a  -  1  -  exp(-2Bk) 

«  ~  exp(2Bh)  -  exp(-2Bk) 


C 


P 


2  '  Bk  Bh  ' 


(11-41) 

(11-42) 

(11-43) 


If  the  node  P  is  placed  on  the  wake  centerline  where 
the  normal  convective  velocity  component  v  is  equal  to 
zero,  then 
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k 

h  +  k 


(11-44) 


A 


S 


h 

h  +  k 


(11-45) 


C 


P 


hk 

2 


(11-46) 


With  the  source  term  linearized  as 

<Vp  ■  ♦p  *  4 

the  following  four  point  partially  parabolic 
discretization  formula  can  be  obtained  by  manipulating 
Eq. (11-40) ,Eq. (11-35)  and  Eq. (11-47) . 


♦p  “  N  *  *s  ‘•’s  *  'Ni  *  ‘’(t. 


(11-48) 


where 


""w  *  S  <Vp 

<  V  p  p 

Ap  =  1  +  Cp  (  (D^)p  +  )  (11-50) 

.n-1  C 

^  “  ^P  ‘  ^  ^  ^ 

and  Aj^,  Ag  and  Cp  are  given  by  Eq.  (11-41)  through 
Eq. (11-46) . 

This  partially  parabolic  form  of  finite  analytic 
formulation  may  be  interpreted  as  a  combination  of  the 
upwind  scheme  in  the  flow  direction  and  the  exponential 


scheme  in  the  normal  direction. 
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II. 5. 2  Elliptic  Form  of  Finite  Analytic  Method  in  the 
Boundary-Fitted  Coordinate  System 


The  non-dimensionized  governing  equations  based  on 
the  cylindrical  coordinate  system  can  be  written  in  the 
transformed  coordinate  system  as  following  using  the 
partial  transformation  where  the  original  orthogonal 
velocity  components  are  left  as  dependent  variables; 


[r(b^u  +  b^v)]^  +  [r(b^u  +  b^v)  =  0 


(11-52) 


9“  -f  2B^4.^  +  (11-53) 


where 

J 


<t>  <|) 


-(j,  -  -4, 

and 


( bi 

u  +  b2  V  )  -  f^ 

(11-54) 

( b^ 

u  +  b2  V  )  -  f^ 

(11-55) 

(11-56) 

^4ti 

(11-57) 

.  1 

b^  =  y, 


Tj.  ^2 


X  b? 

T\f  1 


(11-58) 


J  =  r  (  ^4  Vti  -  ^4  ^ 


(11-59) 
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If  we  locally  linearize  Eq. (11-53)  in  each  numerical 
element,  =  AT|  =  1,  by  the  value  at  calculation  point  P 
as  shown  in  Fig.II-7  and  linearize  the  source  term  by 
Eq. (11-47),  one  obtains  following  equation. 

't’n'l  -  2<Vp  +  G  (11-60) 

where 


^^<|)^P  .n-1  P  .  C 

G  =  (  <l>p  -  <(>p  )  +  <{)^  + 


At 


'<!)  '"p  -(]) 


(11-61) 


If  we  introduce  the  grid  stretching  functions  such 


that 


(11-62) 


then  Eq. (11-60)  can  be  written  in  a  standard  finite 
analytic  form. 


+  4)^*^*  =  2A  4)^*  +  2B  4)^*  +  G 


where 


(11-63) 


(11-64) 


and  the  grid  size  is  stretched  to 
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Equation  (11-63)  is  solved  analytically  using  the 
method  of  separation  of  variable  with  boundary  condition 
expressed  as  the  combinations  of  constant,  linear  and 
exponential  functions  involving  the  nodal  values. 
Evaluation  of  the  analytic  solution  at  the  center  node 
then  provides  a  nine  point  algebraic  discretization 
formula  [3] . 


8 

A  ih 

p 

~  '^nb  ^nb  ^  ^<j> 

n=l 

(11-66) 

where 

^  = 

exp (-Ah)  ^ 

2cosh(Ah) 

(11-67) 

exp  (Ah) 

(11-68) 

2cosh(Ah) 

exp(-B)c)  ^ 

2cosh(Bk)  A 

(11-69) 

exp(Bk) 

(11-70) 

2cosh(Bk)  ^A 

exp(Ah+Bk) 

(11-71) 

^SW  = 

4cosh (Ah) cosh (Bk) 

exp(-Ah+Bk)  ^ 

(11-72) 

""SE  = 

4cosh (Ah) cosh (Bk) 

exp(Ah-Bk)  ^ 

(11-73) 

4cosh (Ah) cosh (Bk) 

^E 

exp(-Ah-Bk)  ^ 

(1-P;,-Pb) 

(11-74) 

4cosh (Ah) cosh (Bk) 

Ap  = 

<^<j)^P  P 

1  ^  > 

(11-75) 
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.n-1  C 


h  tanh(Ah)  ^ 

2A  a' 


(11-76) 


k  tanh(Bk)  ^ 
2B 


(1-Pg)  (11-77) 


and 


=  4Ah  cosh (Ah)  cosh(Bk)  coth(Ah)  E2  (11-78) 


B 


=  1  + 


Bh  coth(Bk) 
Ak  coth(Ah) 


(Pa-1) 


(11-79) 


with 

00  -(-1)"^  h 

E  =  X  - 

"'=l[(Ah)^  +  (X  h)^]^  cosh-^/  (A^  +  B^  +  Xj-  ) y} 
m  T  m 

(11-80) 

X  h  ==  (  m  -  -^  )  7C  (11-81) 

m  z 


For  large  cell  Reynolds  numbers,  the  series  summation 
in  can  be  avoided  by  the  following  asymptotic 
expressions  of  P^  and  P„  based  on  the  theory  of 

A  £3 

characteristics. [4] 


Ak  coth  Ah  >  Bh  coth  Bk; 


P^  =  0,  Pg  =  1  -  (  Bh  coth  Bk  )/(  Ak  coth  Ah  ) • 
Ak  coth  Ah  <  Bh  coth  Bk; 


Pg  =  0,  P^  =  1  -  (  Ak  coth  Ah  )/(  Bh  coth  Bk  )  (11-82) 
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II. 5. 3  Partially  Parabolic  Form  of  Finite  Analytic  Method 
in  the  Boundary-Fitted  Coordinate  System 

The  problems  considered  in  the  present  study  are 
mostly  free  from  separation  and  have  a  predominant  flow 
direction.  In  these  problems,  diffusion  in  the  flow 
direction  is  much  smaller  than  diffusion  in  the  normal 
direction.  Thus  the  diffusion  term  in  the  flow  direction 
can  be  neglected  without  loss  of  accuracy.  This  partially 
parabolic  form  has  advantages  over  the  fully  elliptic  form 
in  the  fact  that  the  computing  time  for  the  calculation  of 
finite  analytic  coefficients  can  be  reduced.  In  the 
situations  when  the  grid  ratio  and  cell  peclet  number 
become  too  large,  the  nine  point  finite  analytic 
coefficients  requires  more  computational  efforts  to  be 
evaluated  accurately.  The  partially  parabolic  form  of 
finite  analytic  method  can  be  used  without  any 
difficulties.  The  partially  parabolic  formulation  can 
also  be  used  for  the  simulation  of  small  separated  flow 
with  FLARE  approximation  [61] . 

Eq. (11-60)  can  be  approximated  using  the  partially 
parabolic  assumption  (  0^^  *  0  )  and  the  upwind  difference 
for  the  convection  term  in  the  ^-direction  as  the 
following  form; 
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22 


<l>Tm  ”  2(B^)p  +  (D^)p  {  <Dp  -  4>w  )  +  G  (li¬ 


es) 


where 


'Vp  ■  "'Vp 


(11-84) 


Introduction  of  the  grid  stretching  functions  T]*  = 


Tl/  \gr  leads  Eq. (11-83)  to  the  standard  finite  analytic 
form. 


22 


*-n*  ~  ^ 


where 


B  = 


<Vp 


V 


22 

gp 


<Vp  ^  ‘''p 


-  Ow  )  +  G 


and  the  grid  size  is  stretched  to 


(11-85) 


(11-86) 

(11-87) 


ATI*  =  k  = 


V 


22 

5p 


(11-88) 


Eqpiation  (11-85)  can  be  solved  analytically  using  the 
north  and  south  nodal  values  as  boundary  conditions. 
Evaluation  of  the  analytic  solution  at  the  center  node 
then  provides  a  four  point  algebraic  discretization 
formula . 


fp  *  ^1  fN  *  '•’s  *  'Hi  ^4 


(11-89) 


where 
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exp (-Bk) 
^  “  2cosh(Bk) 


(11-90) 


exp (Bk) 
2cosh(Bk) 


(11-91) 


b. 

9 


S  <Vp 


1  +  Cp  (  (D^)p  + 


^Vp 

At 


P 

) 


i!$lp  ,n-l 


At 


0 


+  ) 


k  tanh(Bk) 
2B 


(11-92) 

(11-93) 

(11-94) 

(11-95) 


II-6 _ Finite  Volume  Method 


The  strong  conservative  form  of  governing  equations 
in  the  transformed  domain  can  be  written  as  follows; 

(prU)^+(prV)^  =0  (11-96) 

^1  1 

(  p  r  U  0  )^  +  (  p  r  V  (|>  )^  -  t  j  <  ^1  +  ^2 

^2  2 

+  [  <|)^  +  D2  <})^  )  +  J  (11-97) 

where 

U  =  (  bj;  u  +  b2  V  )  (11-98) 

=  (  b^  u  +  b2  V  ) 


V 


(11-99) 
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and 


12  2  2 
Di  =  r  (  ^Tl  ^Tl  ) 

2  2  2  2 
D2  =  r  (  ) 


D 


1 

2 


X-X 

^  T1 


Vtl 


(Il-lOO) 

(II-lOl) 

(11-102) 


In  the  finite  volume  method  [2],  the  governing 
differential  equations  are  integrated  over  a  finite  number 
of  control  volumes  with  size,  =  Atj  =  1,  as  shown  in 
Fig.II-8.  Using  the  Gauss  theorem,  the  volume  integral  of 
terms  under  differential  operator  can  be  converted  into 
surface  integrals  over  the  four  faces  of  a  two  dimensional 
control  volume.  The  resulting  balance  equations  for  each 
control  volume  and  variable  (j>  can  be  expressed  as  follows. 

-  f  J  Sx  d^  dTl  (11-103) 

e  w  n  s  9 

where 

^1  1 

=  (  p  r  U  <|>  )^  -  [  ^  (  <{>^  +  D2  <|)^  )  ]g  (11-104) 
=  (  p  r  V  <>  )^  -  [  ^  (  <|>^  +  (11-105) 


The  volume  integral  of  the  source  term  may  be 
evaluated  as  follows; 


f  J  d£  dh 

Jav  <|>  ^  ' 


P  ^  C 

<  <>P  + 


)  AV 


(11-106) 
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By  putting  the  non-orthogonal  terms  into  source  term, 
the  Eq. (11-103)  can  be  written  as  as  follows. 

r  r 


^  2 

+  [prV<|)-^D2<|>^]^-[prV 


2 

T  4  ^ 


(11-107) 


where  is  the  source  term  arised  from  the  non- 
orthogonality  of  numerical  grid  and  defined  as; 

b^Pi  ^1  ^2  2 

%  '  '  J  ‘>2  J  “2  fTl>w+<  J  ‘’l  Vn-<  J 


(11-108) 


The  Eq. (11-107)  is  the  standard  finite  volume 
formulation  and  discretization  of  this  equation  is  well 
explained  in  Patankar  [2] .  Using  the  hybrid  numerical 
scheme,  the  discretization  equation  can  be  obtained  as 
follows . 


*p  ♦p  '  *E  ♦w  N)  ♦n  *S  *3  “’ll) 


where 


Ag  =  MAX  [ 

A^^  =  max  [ 

Ajj  =  MAX  [ 


0.5*1  I  , 

0.5*1  I  , 

0.5*1  F^  I 


D  ]  -  0.5*F^ 
e  e 

D  ]  +  0.5*F 
w  w 

Dn  ]  -  0.5*F^ 


(11-109) 

(II-llO) 

(II-lll) 


r 


(11-112) 
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= 

MAX  [  0 

.5*1  Fg  1  ,  Dg  ]  +  0.5*Fg 

(11-113) 

(11-114) 

c  ^ 

Sx  Av  + 

b 

(11-115) 

and 

F  = 
e 

(  p  r  U 

)e 

(11-116) 

F  = 
w 

(  p  r  U 

(11-117) 

F  = 
n 

(  p  r  V 

)n 

(11-118) 

^s  “ 

(  p  r  V 

)s 

(11-119) 

with 

°e  = 

1 

<  T°i 

^e 

(11-120) 

1 

<  j  °i 

)w 

(11-121) 

D  = 
n 

2 

<  J  ^2 

>n 

(11-122) 

D  = 
s 

2 

<  J  °2 

^s 

(11-123) 

In  these  equations,  the  notation  MAX  [  A  ,  B  ]  means 
the  greater  of  A  and  B  and  |  A  |  means  the  absolute  value 
of  A. 


The  problems  considered  in  the  present  study  are 
mostly  free  of  separation.  Thus  the  numerical  false 
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diffusion  problem  has  a  negligible  effect  on  the  accuracy 
of  solution.  In  these  situations,  the  choice  of  numerical 
scheme  is  not  important.  However,  the  use  of  different 
numerical  method  requires  different  usage  of  wall  function 
method  in  the  turbulent  flow  calculation.  The  finite 
analytic  method  usually  uses  the  two-point  wall  function 
method  while  the  wall  function  method  of  Launder  and 
Spalding  [43]  is  employed  for  most  of  finite  volume 
calculations.  The  primary  numerical  scheme  for  the 
present  study  is  the  finite  analytic  method.  However,  the 
finite  volume  method  is  employed  when  the  use  of  the  two- 
point  wall  function  method  is  restricted  by  a  certain 
problem  or  when  the  comparisions  of  predictions  by  two 
different  methods  are  thought  to  be  important.  It  is 
noted  that  the  finite  analytic  method  usually  employs  non- 
dimensionized  form  of  governing  equations  while  finite 
volume  method  employes  dimensional  form.  All  the 
derivations  follow  this  common  practice.  The  use  of  each 
numerical  scheme  for  a  certain  problem  is  explained  when  a 
specific  problem  is  mentioned. 
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CHAPTER  III 
SOLUTION  PROCEDURE 

III.l  Introduction 

A  difficulty  encountered  in  solving  the  Navier-Stokes 
equations  for  incompressible  flow  concern  the  handling  of 
pressure  terms.  The  difficulty  arises  from  the  fact  that 
pressure  does  not  have  its  own  equation  in  the  sets  of 
Navier-Stokes  equations.  The  pressure  term  only  appears 
in  the  momentum  equations  which  must  be  considered  as  the 
governing  equations  for  velocity  components.  Thus  the 
continuity  equation  should  be  used  to  solve  the  pressure 
variable.  However,  there  is  no  pressure  term  in  the 
continuity  equation.  This  difficulty  is  usually 
circumvented  by  deriving  the  pressure  equation  by 
manipulating  the  continuity  equation  and  momentum 
equations.  It  is  noted  that  the  algebraic  relations 
obtained  in  the  previous  chapter  can  be  used  only  for  the 
velocity  components  and  turbulent  quantities.  The  way  of 
resolving  pressure  varible  and  satisfying  mass 
conservation  should  be  sought.  In  the  present  study,  the 
SIMPLER  (  Semi-Implicit  Method  for  Pressure  Linked 
Equation  -  Revised  )  numerical  algorithm  of  Patankar  [2] 
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is  imployed  for  this  purpose.  In  this  algorithm,  a 
Poisson  like  equation  which  is  derived  from  the 
manipulation  of  continuity  equation  and  momentum  ecjuations 
is  solved  to  obtain  the  pressure  field  and  in  addition, 
the  pressure  correction  equation  is  again  solved  to 
correct  velocity  components  to  satisfy  continuity 
equation.  Following  sections  are  devoted  to  the 
derivation  of  pressure  and  pressure  correction  equation  in 
the  Cartesian  coordinate  system  and  in  the  body-fitted 
coordinate  system. 


III. 2.  Pressure  and-Prassure  Correction .Equation  in 

Cartesian  Coordinate  System 


From  the  discretization  equations  for  momentum 
equations  such  as  Eq. (11-40)  or  Eq. (11-66),  the  velocity 
field  (  u  ,  V  )  can  be  decomposed  into  pseudovelocity 
components  (  u  ,  v  )  and  the  pressure  gradient  terms 
contained  in  the  source  term  as; 


u  = 

A 

u 

( 

Pc 

~  P  ) 

e 

e 

e 

P 

V  = 

A 

V 

d  ( 

P»T 

-  P  ) 

n 

n 

n 

P 

where 

A 

A 

1 

^e 

1 

<  I 

nb 

'  I 

nb 

A  V. 

nb 

u  .  +  b 

nb  u 

V  = 

n 

A 

n 

A  * 
nb 

V  .  +  b 

nb  V 

(III-l) 

(III-2) 

)g  (III-3) 


(III-4) 
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and 

■le  =  '  ‘«u'e  W  [  (5x)^  1  (III-5) 

-In  -  '  ‘"v'n  >  ^  ^ 

In  these  equations,  subscripts  e  and  n  denote  the 
value  evaluated  at  staggered  velocity  nodes  e  and  n  as 
shown  in  Fig.III-1.  For  example  A^,  are  algebraic 
coefficients  Ap,  Cp  evaluated  at  node  e  and  (Sx)^,  (Sy)^ 

are  the  distances  between  pressure  calculation  points. 

The  variables  R^,  R^  are  the  effective  Reynolds  number  for 
variable  u  and  v.  All  the  coefficients  are  based  on  the 
finite  analytic  method  in  non-dimensional  form  given  in 
Eq. (11-41)  through  Eq. (11-43) . 

The  continuity  equation  within  the  control  volume 
cell  can  be  wriiten  as  following; 

(Ay)eUg  -  (Ay)^u^  +  (Ax)jjV^  -  (Ax)  ^v^  =  0  (III-7) 

where  Ax,  Ay  are  the  sizes  of  control  volume  cell  as 
shown  in  Fig.III-1. 

The  pressure  equation  can  be  derived  by  substituting 
momentum  equation,  Eq. (III-l)  and  Eq. (III-2)  into 
continuity  equation,  Eq. (III-7) . 

^  ^  ^E  ^  ^  ^  ^  ~  ^  (III-8) 
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where 

Ae  =  e 

(III-9) 

A„  -  (Ay)„  d„ 

(III-IO) 

\  ■  <'^='>1,  ‘‘n 

(III-ll) 

As  =  (Ax)^ 

(III-12) 

(III-13) 

and 

D  =  (Ay)^u^  -  (Ay)^u„  +  (Ax)^v^  -  s^s 


The  velocity  components  (  u  ,  v  )  obtained  from  the 
solution  of  the  momentum  equations  will  generally  not 
satisfy  the  continuity  equation  unless  the  pressure  field 
is  correct .  We  denote  these  imperfect  velocity  components 
as  starred  velocity  components  (  u*  ,  v*  ) .  These  starred 
velociy  components  must  be  corrected  to  satisfy  continuity 
equation . 


u 


e 


r  f 


(III-15) 


V 

n 


n 


(III-16) 


t 


where  P  is  the  pressure  correction. 
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The  pressure  correction  equation  is  obtained  by 
substituting  these  velocity  correction  equations  into 
continuity  equation,  Eq. (III-7) . 

Pp  -  'S;  4  *  4  *  Si  '■n  -  (Ill-n) 

where 

D*  =  (Ay)^Ug  -  (Ay)^u„  +  (Ax)^Vj^  -  (Ax)  (III-18) 

and  the  coefficients  A^,  A^,  A^,  Ag  and  Ap  are  same 

★ 

as  those  given  in  pressure  equation  and  the  mass  source  D 
is  based  on  starred  velocity  components  which  are  obtained 
from  the  solution  of  momentum  equations. 

After  the  pressure  correction  equation  is  solved,  the 
velocity  components  are  corrected  through  velocity 
correction  equations  to  satisfy  continuity  equation. 

III. 3  Pressure  and  Pressure  Correction  Equation  in 
Boundary-Fitted  Coordinate  System 

III. 3.1  One  Velocity  Staggered  Grid  Method 

In  one  velocity  staggered  grid  method,  only  one 
velocity  component  is  stored  in  each  staggered  node  e  and 
n  as  shown  in  Fig.III-2.  Thus  derivation  of  pressure  and 
pressure  correction  equation  in  this  grid  arrangement  is 
not  much  depart  from  those  in  Cartesian  coordinate  system. 
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The  discretized  monentum  equations  for  and  can 
be  written  as  follows; 


u  = 
e 

(H 

u 

>e  - 

E-Pp 

)  - 

(D 

2)e(P 

-p  ) 
en  ^es 

(III-19) 

V  = 

n 

(H 

V 

>n  - 

(D^)n(P 

“P 
en  ^ 

wn^ 

- 

^^I)n 

<PN-Pp) 

(III-20) 

where 

(H  ) 
u  e 

=  ' 

e 

I  V  “nb 
nb 

)e 

(III-21) 

(H  ) 

V  n 

=  ' 

i  < 

n 

I  ''nb  ^ 

nb 

b 

V 

)n 

(III-22) 

and 

<°l'e 

= 

f  ^e 

b\  ) 

e  ^ 

/ 

(III-23) 

<“2>e 

= 

t  ^e 

b^ ) 

e  1 

/ 

^e  ^ 

(III-24) 

= 

^  ^n 

(  R  r 

V 

b^  ) 

n  ^ 

/ 

^n  1 

(III-25) 

= 

^  ^n 

(  R  r 

V 

bl  ) 

n  1 

/ 

^n  ] 

(III-26) 

In  these 

equations. 

the 

algebraic  coefficients  are 

.  on  the 

finite  analytic 

method 

with 

non-dimensional 

given 

in 

chapter  II- 

■4. 

The  velocity  components  (  u  ,  v  )  are  then  decomposed 

A  A 

into  pseudovelocity  components  (  u  ,  v  )  and  the  pressure 
gradient  terms . 
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(III-27) 


"n  =  "n  -  (°2)n  <  ^N  “  ^P  > 


(III-28) 


where 


u  =  (H  )  -  (D,)  (  p  -  p  ) 

e  u  e  '  2'e  ^en  ^es 


(III-29) 


V  =  (H  )  -  (D,  )  (  p  -  p  ) 

n  V  n  In  ^en  ^wn 


(111-30) 


The  cross  derivative  terms  of  pressure  are  contained 
in  the  pseudovelocity  components  to  avoid  a  complicate 
nine  point  pressure  equation  and  are  explicitly  evaluated 
using  the  pressure  obtained  from  the  previous  iteration. 


The  continuity  equation  within  the  control  volume 
cell  in  non-orthogonal  coordinate  system  can  be  written  as 
follows . 

[r(bju  +  b2v)]g  -  [r(bj^u  +  bjv)]^ 

+  [r(b^u  +  b^v)]^  -  [r(b^u  +  b^v) ]  =  0  (III-31) 


The  Eq. (III-31)  can  be  rearranged  as  following  form. 


(rbiu)e  -  {rh^u)^  +  (rb2v)^  -  (rb2v)g  = 


-  D„  (III-32) 


where 


Dn  =  (rb2V)^ 


(rb2V)w  +  (rb^u)^  -  (rb^u)^ 


(III-33) 
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The  pressure  equation  can  be  derived  by  substituting 
Eq. (III-27)  and  Eq. (III-28)  to  continuity  equation, 

Eq. (111-32) . 


^  Pp  -  'S;  Pe 

(III-S4) 

*E  ■  <  "  Pi  >e 

K'e 

(III-S5) 

Pi  >w 

(III-S6) 

'Si  -  <  P2  'n 

(III-S7) 

"  P2  <3 

<°2'3 

(III-S8) 

(III-S9) 

where 

D  =  (rb];u)^  -  (rbj^u)^  +  (rb^v)^  -  (rb^v)^  + 

(III-40) 

A 

and  the  mass  source  D  is  based  on  the  pseudovelocity 
components.  The  mass  source  term  Djj  expressed  in  Eq.  (HI¬ 
SS)  should  be  evaluated  from  the  continuity-satisfying 
velocity  components  from  the  previous  iteration  level.  It 
is  noted  that  the  velocity  components  that  are  needed  for 
evaluation  of  Dj^  are  not  stored  at  the  staggered  nodes  and 
are  obtained  by  averaging  neighbouring  nodal  values . 


The  imperfect  velocity  components  {  u  ,  v  )  which 
are  obtained  from  the  solution  of  momentum  equations  are 
corrected  by  the  pressure  correction  in  order  to  satisfy 
the  continuity  equation. 

“e  -  “e  -  -  f'p  >  (III-41) 

''n  -  ''n  -  <°2>n  <  '  '’p  >  (III-42) 

I 

where  P  is  the  pressure  correction.  The  pressure 
correction  terms  arised  from  the  cross  derivatives  are 
neglected  to  ensure  the  diagonal  dominant  pressure 
correction  equation  which  is  an  important  parameter  to 
achieve  the  stability  of  solution  procedure.  The 
neglection  of  these  terms  do  not  effect  the  accuracy  of 
final  converged  solution  since  the  pressure  correction 
should  be  zero  at  the  final  converged  stage.  However, 
this  practice  may  slow  down  the  convergence  when  the 
numerical  grid  is  strongly  non-orthogonal . 

The  pressure  correction  equation  is  obtained  by 
substituting  velocity  correction  equations  into  continuity 
equation,  Eq. (III-32) . 

*P  ^P  •  »E  "e  ^  fw  *  "“n  *  *^3  (111-13) 


where 
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D  =  (rb]^u*)g  -  (rb^u*)^  +  (rb^v*)^  -  (rb^v*)^  +  Dj^ 

(III-44) 

where  the  coefficients 
same  as  those  given  in  pressure  equation  and  D*  is  based 
on  the  starred  velocities  {  u*,  v*  )  which  only  satisfy 
the  monentum  equation. 

With  the  updated  pressure  correction,  the  velocities 
are  improved  through  the  velocity  correction  equations  to 
satisfy  continuity  equation. 

It  is  noted  that  pressure  and  velocity  coupling 
technique  used  in  the  present  study  is  slightly  different 
from  that  used  by  Shyy  et  al.  [24]  although  same  one 
velocity  staggered  grid  method  is  employed. 

III. 3. 2  Two  Velocities  Staggered  Grid  Method 

In  two  velocities  staggered  method,  both  velocity 
components  are  stored  at  each  staggered  node  e  and  n  as 
shown  in  Fig.III-3.  The  discretized  form  of  momentum 
equations  for  these  velocity  components  can  be  written  as 
follows  with  the  pressure  gradient  term  expressed 
explicitly. 

“e  -  '•’uU-  <°l>e  <PE-Pp>  -  <P2>e  'Pen-Pes' =  > 


'  *N'  *S 
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=  <«v)e  -  "  ^^2^  ^Pen'Pes^ 
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(III-51) 
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(III-52) 
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(III-53) 
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(III-54) 

In  these  equations,  the  algebraic  coefficients  are 
based  on  finite  volume  method. 


The  velocity  field  (  u  ,  v  )  can  be  decomposed  into 

A  A 

pseudovelocity  components  (  u  ,  v  )  and  the  pressure 
gradient  terms  as  follows. 


(D^)e  (  Pe  -  '’P  ) 


(III-55) 
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(III-57) 
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(III-61) 

(III-62) 


The  cross  derivative  terms  of  pressure  are  again 
contained  in  pseudovelocity  components  to  obtain  diagonal 
dominant  five  point  pressure  ecjuations. 


The  continuity  equation  for  the  control  volume  cell 
shown  in  Fig.III-3  can  be  written  as  follows; 

U  -U  +V  -V  =0  (III-63) 

e  w  n  s 

where  U  and  V  are  the  contravariant  velocity 
components  defined  as; 
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U  =  r  {  u  +  b2  V  ) 

(111-64) 

2  2 

V  =  r  (  bj^  u  +  b2  V  ) 

(III-65) 

The  pressure  equation  is  then  obtained  by 

sustituting 

Eq. (111-55)  through  Eq. (III-58)  into  continuity 

Eq. (III-63) . 

equation, 

Ap  Pp  -  Ae  Pe  \  Pm  'Si  Pn  Ps  '  “ 

where 

(III-66) 

Ae  ’  Pe  (  d"  + 

(III-67) 

'H.  -  r„  (  bi  d“  +  hi  dJ  )„ 

(III-68) 

-  r„  (  hi  dP  *  hi  dJ  )„ 

,, 2  u  ^2  v 

(III-69) 

*3  “  '  bp  ^2  b2  ^2  '  s 

(III-70) 

A  A  A  A  A 

(III-71) 

D  =  U  -  U  +  V  -  V 
e  w  n  s 

A  A 

(III-72) 

and  U  and  V  are  based  on  the  pseudovelocity 

components  and  defined  as; 

A  1  1 

U  =  r  (  b^  u  +  b2  V  ) 

A  2  ^  2  ^ 

V  =  r  (  b^  u  +  b2  V  ) 

(III-73) 

(III-74) 

(III-74) 
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Since  the  velocity  components  obtained  from  the 
solution  of  momentum  equations  with  imperfect  pressure 
field  generally  do  not  satisfy  continuity  ec[uation.  These 
starred  velocity  components  (  u*,  v*  )  are  corrected 
through  the  pressure  correction  in  order  to  satisfy  mass 
conservation  during  the  iteration  process . 


u  = 
e 

★ 

^e  - 

Pp)  - 

(D^)e 

(P 

1  1 

-P  ) 
en  es' 

(III-75) 

V  = 

e 

* 

^e  - 

K>e 

Pp)  - 

<'>2>e 

(P 

1  1 

-P  ) 
en  es' 

(III-76) 

and 

'^n  = 

* 

'^n  - 

^^en 

-Pwn> 

-  (Dj) 

n 

(PN-Pp) 

(III-77) 

II 

c 

> 

* 

^n  - 

<'n 

^^en 

-^wn) 

-  (Dj) 

n 

(Pn’^p) 

(III-78) 

where  P  is  the  pressure  correction 

In  order  to  facilitate  the  velocity  correction 
procedure,  following  contravariant  velocity  correction 
equations  are  derived  using  Eq. (III-75) -Eq. (III-78)  and 
Eq. (III-64)-Eq. (111-65) . 


U 

e 

it 

=  Ue  - 

(  P£  ~  Pp  ) 

(III-79) 

V 

n 

it 

=  Vn  - 

(III-80) 

where 
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*  1  *  1  * 
Ue  =  Te  (  bi  u  +  b2  V 


''n  -  *  *>2  'n 


and 


,  1 

u 

,  1 

'tt)  =  ^ 

U'  e  e 

( 

^1 

^1 

+ 

^2 

2 

u 

2 

=  r 

V'  n  n 

( 

°2 

+ 

b2 

(III-81) 

(III-82) 

(III-83) 

(III-84) 


The  pressure  correction  terms  arised  from  the  non¬ 
orthogonality  of  numerical  grids  are  neglected  to  avoid 
the  complicated  nine  point  pressure  correction  equation 
and  to  ensure  diagonal  dominant  pressure  correction 
equation . 


The  pressure  correction  equation  is  obtained  by 
substituting  these  contravariant  velocity  correction 
equations  into  continuity  equation,  Eq. (III-63) . 

'■p  -  *E  Pe  *  '■w  +  '■n  *  '‘S 

where 

-  Oe  -  “w  +  V*  -  V*  (III-86) 


where  the  coefficients  A^^,  Ag  and  Ap  are  the 
same  as  those  in  pressure  equation  and  D*  is  based  on  the 
starred  velocities  which  are  obtained  from  the  solutions 


of  the  momentum  equations. 
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From  the  updated  pressure  correction,  the 
contravariant  velocity  components  are  corrected  to  satisfy 
mass  conservation.  The  other  contravariant  velocity 
components,  and  which  lie  parallel  to  the  control 
volume  surfaces  are  obtained  from  linear  interpolation  of 
newly  computed  neighbouring  nodal  values . 

Then  the  corrected  velocities  are  obtained  from  these 
continuity-satisfying,  corrected  contravariant  velocity 
components . 

u=  (b2U-b2V)/J  (III-87) 

v  =  -  (b^U-bj;v)  /  J  (III-88) 

where 

J  =  r  (  h\  bl  -  hi  hi  )  (III-89) 

It  is  noted  that  the  conitinuity-satisfying, 
corrected  contravariant  velocity  components  are  stored  and 
used  for  the  evaluation  of  algebraic  coefficients  A^,  A^, 
A^,  Ag  and  Ap  for  transport  quantities. 

Maliska  and  Raithby  [30]  used  the  same  procedure  as 
outlined  above  in  the  calculation  of  three  dimensional 
parabolic  flow  of  arbitrary  cross  sectioned  ducts. 

However,  they  used  exact  nine-point  pressure  and  pressure 
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correction  equations  which  may  deteriorate  the  stabillity 
of  solution  process  when  the  numerical  grids  are  strongly 
non-orthogonal . 

Ill .4  Overall  Solution  Procedure 

The  solution  procedure  of  SIMPLER  algorithm  is 
adopted  for  all  the  calculations  in  the  present  study. 

The  momentum  equations  and  turbulent  transport  equations 
are  first  solved  by  parabolic  marching  technique  in  the 
flow  direction  using  the  pressure  obtained  from  the 
previous  iteration.  After  che  marching  of  calculations 
for  turbulent  quantities  is  completed,  the  pressure  field 
is  solved  elliptically  with  several  global  iterations. 

The  advantage  of  this  global,  elliptic  pressure 
calculation  technique  is  a  monotonic  and  rapid  convergence 
of  pressure  due  to  a  proper  consideration  of  the  elliptic 
nature  of  the  pressure  field,  especially  near  the  leading 
and  trailing  edges . 

The  detailed  overall  solution  procedure  can  be 
outlined  as  following. 

1)  Generate  the  numerical  grids  using  the  methods 
outlined  in  chapter  II-3  and  save  the  geometric  data,  such 
as  f  /  9  etc. 


2)  Specify  the  inlet  boundary  conditions  for 
transport  cjuantities,  the  velocity  components  and 
turbulent  quantities. 

3)  Set  pressure  and  pressure  correction  equal  to  zero 
everywhere  in  the  solution  domain. 

4)  Calculate  the  algebraic  coefficients  and  source 
terms  for  momentum  equations,  such  as  those  given  in 

Eq. (11-67)  through  Eq. (11-81) ,  at  downstream  station  using 
the  transport  quantities  obtained  from  the  previous 
iteration  level.  This  step  includes  the  calculation  of 
coefficients  of  pressure  and  pressure  correction  equation 
such  as  those  given  in  Eq. (III-35)  through  Eq.  (III-39)  . 

5)  Solve  the  momentum  equations,  such  as  Eq. (I 11-19) 
and  Eq.(III-20),  using  the  pressure  from  the  previous 

A  if 

iteration  level  to  obtain  starred  velocity  field,  (  u  ,  v 

)  . 

6)  Using  this  starred  velocity  components,  calculate 

* 

the  mass  source  for  the  pressure  correction  equation,  D  . 

7)  Solve  the  pressure  correction  equation,  such  as 
Eq. (III-43)  and  update  the  velocity  field  using  Eq. (III- 
41)  and  Eq. (III-42)  to  satisfy  continuity  equation. 
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8)  Calculate  the  pseudovelocity  field,  (  u  ,  v  )  from 
Eq. (III-29)  and  Eq. (III-30)  using  the  corrected  velocity 
components . 

9)  Calculate  the  mass  source  for  pressure  equation, 

yv 

D,  and  store  for  later  use. 

10)  Calculate  the  algebraic  coefficients  and  source 
terms  for  the  turbulent  transport  equations  using  the 
updated  velocity  field  and  turbulent  quantities  from  the 
previous  iteration  level . 

11)  Solve  the  turbulent  transport  equations. 

12)  March  to  the  next  downstream  station  and  repeat 
the  procedure  4)  to  11) . 

13)  After  the  marching  is  completed  at  the  last 
downstream  station,  the  pressure  equation  such  as  Eq. (III- 
34)  is  solved  elliptically  with  several  global  iterations 

(  typically  20  sweeps  )  from  downstream  to  upstream  to 
update  new  pressure  field. 

14)  Return  to  step  4)  for  the  next  iteration  level. 

15)  Steps  4)  through  14)  are  repeated  until 
convergence  or  steady  state  solution  is  reached. 
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All  the  algebraic  equations  are  solved  by  line  by 
line  technique  using  tridiagonal  matrix  algorithm.  It  is 
noted  that  although  parabolic  marching  technique  is 
employed  for  the  calculation  of  transport  quantities,  all 
the  transport  quantities  are  stored  in  two  dimensional 
array.  Thus  the  present  solution  procedure  has  potential 
for  the  calculation  of  separated  flow. 
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CHAPTER  IV 

CALCULATION  OF  TURBULENT  WAKE  PAST  A  FLAT  PLATE  BY  WAKE 

FUNCTION  METHOD 

IV.  3 _ Introduction 

The  turbulent  flow  past  a  finite  flat  plate  has  been 
studied  by  many  researchers  because  it  provides  a 
fundamental  understanding  of  wake  development  and  the 
basic  feature  of  viscous-inviscid  interaction  at  the 
trailing  edge  of  the  body.  According  to  Alber  [62],  the 
wake  region  can  be  divided  into  three  fundamental  regions; 
laminar  wake  region  where  the  laminar  sublayer  on  the 
plate  is  destroyed,  turbulent  inner  near  wake  region  where 
the  logarithmic  remnant  of  the  trailing  edge  boundary 
layer  is  destroyed,  the  far  wake  region  where  the  flow 
field  loses  the  memory  of  the  turbulent  boundary  layer  on 
the  flat  plate  and  attains  the  self  preserving  form.  The 
experimental  data  of  Chevray  and  Kovasznay  [63],  Ramaprian 
et  al.  [64],  Pot  [65],  Andreopoulous  and  Bradshaw  [66] 
confirm  such  a  wake  behaviour  behind  a  flat  plate.  Some 
comparisions  and  reviews  of  these  data  can  be  seen  in 
Ramaprian  et  al.  [64]. 
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The  earlier  calculations  by  Rodi  [67],  Launder  et 
al.  [68]  had  placed  more  emphasis  on  the  test  of  turbulence 
model.  More  extensive  calculations  were  carried  out  by 
Patel  and  Scheuerer  [69]  for  symmetric  and  asymmetric 
turbulent  wake  of  a  flat  plate.  The  results  were  compared 
with  several  available  experimental  data.  In  these  works, 
the  authors  used  the  boundary  layer  calculation  method 
with  zero  pressure  gradient  and  calculations  were  started 
at  the  trailing  edge  of  the  flat  plate  with  initial 
conditions  provided  by  experimental  data.  Thus,  the 
boundary  layer  on  the  flat  plate  and  viscous-inviscid 
interaction  at  the  trailing  edge  of  flat  plate  were  not 
considered.  A  detailed  large  domain  solution  surrounding 
the  trailing  edge  of  the  flat  plate  was  obtained  recently 
by  Patel  and  Chen  [70] .  The  two-layer  model  was  adopted 
for  the  boundary  layer  calculation  on  the  flat  plate  and 
the  k  -  e  turbulence  model  was  adopted  for  the  wake 
calculation.  The  boundary  layer  on  the  flat  plate  was 
accurately  calculated  to  the  lar.in.^r  sublayer  (  y"*"  ~  0.4  ) 
by  the  two-layer  model  in  which  the  turbulent  kinetic 
energy  and  turbulent  kinetic  energy  dissipation  rate  in 
the  laminar  sublayer  and  the  buffer  layer  and  a  part  of 
logarithmic  layer  are  specified  by  universal  functions 
based  on  the  experimental  data  curve-fitting  and  eddy 
viscosity  in  these  regions  is  calculated  by  the  Van  Driest 
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formula.  Comparisons  of  calculated  results  with 
experimental  data  and  results  from  the  previous 
calculations  were  made.  However,  a  detailed  destruction 
of  laminar  sublayer  and  buffer  layer  in  the  wake  was  not 
calculated  since  the  k  -  e  turbulence  model  was  adopted  in 
the  wake  calculation.  For  example,  the  destruction  of 
turbulent  kinetic  energy  peak  at  y"*"  ~  15  on  the  trailing 
edge  of  turbulent  boundary  layer  and  the  development  of 
wake  centerline  velocity  in  the  laminar  wake  region  were 
not  calculated.  While  the  two  layer  model  predicted  a 
reasonably  accurate  result,  it  requires  many  computational 
nodes  in  the  near  wall  region  over  the  plate.  It  is  quite 
questionable  that  the  present  low  Reynolds  number 
turbulence  model  can  be  used  confidently  for  the 
prediction  of  rapid  mixing  in  the  laminar  wake  region 
although  it  has  been  successful  in  the  various 
calculations  of  boundary  layer  type  flow. 

The  purpose  of  the  present  study  is  to  develop  a 
calculation  method  in  which  the  concept  of  a  wall  function 
method  is  adopted  and  extended  to  the  wake  calculation. 

In  this  approach,  the  laminar  wake  region  is  excluded  from 
the  calculation  domain  and  wake  functions  are  introduced 
to  accurately  account  for  the  large  streamwise  velocity 
gradient  along  the  wake  centerline  in  the  same  spirit  that 
wall  functions  are  used  in  the  near  wall  calculation.  It 
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is  shown  that  the  application  of  the  combination  of  wall 
and  wake  function  may  satisfactorily  predict  the  turbulent 
boundary  layer  on  the  plate,  the  trailing  edge  interaction 
and  the  development  of  turbulent  wake. 

IV .2  Numerical  Grids  and  Solution  Domains 

The  solution  domain  shown  in  Fig.IV-1  consists  of 
wall  function  region,  wake  function  region,  and 
calculation  region.  The  wall  function  region  is  one 
computational  control  volume  over  the  plate  which  extends 
from  y"*"  =  0  to  y"*"  '*150.  The  wake  function  region  is  one 
computational  control  volume  just  downstream  of  the 
trailing  edge  of  flat  plate  extending  approximately  from 
x"*"  =  0  to  x'^  ~  500  in  the  wake  region . 

Non-uniform  57x31  grids  are  generated  numerically  by 
the  grid  generation  technique  outlined  in  chapter  II-4-1 
with  exponential  distribution  along  the  y-direction  and 
sinusoidal  distribution  along  the  x-direction  so  that  an 
appropriate  grid  concentration  can  be  made  close  to  the 
wall  and  near  the  trailing  edge  of  flat  plats  as  shown  in 
Fig.II-3.  The  details  of  informations  on  grid  and  sizes 
of  solution  domains  are  given  in  Table  2.  The  numerically 
generated  coordinate  lines  are  treated  as  control  volume 
lines  and  grid  lines  are  placed  at  the  center  of  control 
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volume  lines  as  suggested  by  Patankar  [2] .  Of  the  31  grid 
points  in  the  y-direction,  10  grid  points  are  placed 
within  the  boundary  layer  on  the  body.  The  origin  of 
coordinate  system  is  placed  at  the  trailing  edge  of  plate. 


IV. 3 _ Boundary.  Conditions 


The  upstream  inlet  condition  for  u-velocity  is 
specified  by  the  following  simple  flat-plate  correlations 
[71]  . 

u  =  (  f  (IV-1) 

o 

-1/5 

8(x)  «  0.37  X  (  Re  X  )  (IV-2) 


where  8  (x)  denotes  the  boundary  layer  thickness  at 
distance  x  from  the  leading  edge  and  Re  is  the  Reynolds 
number  based  on  the  free  stream  velocity  and  plate  length. 


The  upstream  inlet  condition  for  turbulent  kinetic 
energy  and  the  rate  of  turbulent  kinetic  energy 
dissipation  are  specified  by  following  simple  relations. 


k  =  c 


-1/2 


e  =  c 


3/4  ,  3/2 


“x  <  ^  8  > 


/  1 


m 


(IV-3) 


(IV-4) 
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where  is  the  friction  velocity,  c^  = 

0.09  and  the  length  scale  Im  is  calculated  by  Escudier 
formula  [72] . 


y/5  <  c^/K  ,  im  -  K  y 
y/5  >  c^/K  ,  Im  =  0^8 


In  the  wall  function  region,  u,  k,  and  E  are 
calculated  by  the  following  conventional  relations 


+ 

y 


+  B 


e  =  u^  /  Ky 

where  y^  *  Re  u  y,  K  is  von-Karman  constant, 
X 


and  B  *  5.5. 


(IV-5) 

(IV-6) 

[43]  . 
(IV-7) 
(IV-8) 
(IV-9) 

K  -  0.42 


The  two  point  wall  function  method  used  in  the 
present  study  is  a  slightly  different  from  the  wall 
function  method  generally  used  in  the  control  volume 
method.  It  is  assumed  that  at  least  two  u  calculation 
nodes  are  placed  in  the  logarithmic  region.  The  friction 
velocity  u^  is  first  determined  by  iterative  method  using 
the  u  value  in  the  third  node  (  y"*"  „  150  )  and  Eq.  (IV-7)  . 
Then  the  u,  k,  e  values  in  the  second  node  (  y"^  60  )  are 


evaluated  using  this  friction  velocity  and  Eqs . (IV-7) - (IV- 
9)  . 


In  the  wake  function  region,  the  wake  centerline 
velocity  u^,  and  other  variables,  v^,  k^,  and  are 
specified  at  the  locations  as  shown  in  the  Fig.IV-2. 
Detailed  specification  and  explanation  of  these  wake 
functions  will  be  given  later. 

The  other  boundary  conditions  are  sama  as  those 
reported  in  chapter  II-3-2. 

iy.4..-Wake  Function  Method 

IV. 4.1.  Significance  of  Wake  Function  Method. 

One  of  the  difficulties  that  arises  in  the  prediction 
of  wake  flow  is  that  when  the  turbulent  boundary  layer 
leaves  the  trailing  edge  of  flat  plate,  the  turbulent 
structure  of  flow  changes  rapidly  because  the  no  slip 
condition  on  the  plate  abruptly  changes  to  the  symmetry 
condition  in  the  wake.  There  is  a  rapid  mixing  near  the 
wake  centerline  and  the  turbulent  boundary  layer  on  the 
plate  breaks  down  and  developes  into  the  wake.  The  flow 
development  in  the  wake  can  be  divided  into  three 
fundamental  regions  as  mentioned  before.  The  practical 
difficulty  encountered  in  the  calculation  of  turbulent 
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wake  is  that  none  of  the  present  turbulence  models  can  be 
used  confidently  for  the  calculation  of  the  laminar  wake 
and  far  wake .  On  the  other  hand,  it  is  not  easy  to 
capture  numerically  the  rapid  variation  of  turbulent 
structure  near  the  wake  centerline  unless  large  number  of 
nodes  are  placed  near  the  trailing  edge.  The  same 
difficulties  are  also  encountered  in  the  calculation  of 
flow  near  the  wall.  The  wall  function  is  introduced  to 
circumvent  these  two  difficulties.  In  the  present  study, 
we  introduce  wake  functions  which  may  account  for  the 
development  of  velocity  along  the  wake  centerline  while 
laminar  wake  region  is  excluded  from  the  calculation 
domain.  The  detailed  calculation  in  the  laminar  wake 
region  and  the  uncertainty  of  the  turbulence  model  in  the 
laminar  wake  calculation  can  be  removed  by  this  approach. 
If  there  is  no  separation  near  the  trailing  edge  of  the 
body  like  present  problem,  the  downstream  near  wake 
calculation  is  strongly  depends  on  the  flow  condition 
immediately  downstream  of  the  trailing  edge.  Thus  a 
proper  imposition  of  the  wake  function  may  reduce  the 
error  initiated  by  the  wake  inlet  condition. 

One  of  the  difficulties  in  the  calculation  of  a 
symmetric  wake  when  the  wall  function  method  is  used  in 
the  calculation  of  boundary  layer  of  body  is  the 
imposition  of  symmetry  condition  for  the  u-velocity  along 
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the  wake  centerline  in  the  turbulent  inner  near  wake 
region.  In  the  initial  phase  of  this  study,  an 
interpolation  function  which  satisfies  the  symmetry 
condition  along  the  wake  centerline  was  used  for  the 
evaluation  of  quantities  of  the  first  control  volume  in 
the  wake  and  along  wake  centerline  while  calculations  were 
carried  out  from  the  second  control  volume .  The  square 
interpolation  function,  Ay2  +  B,  was  used  for  the 
evaluation  of  u,  k  ,  e  and  linear  interpolation  function. 
Ay,  was  used  for  the  evaluation  of  v.  However,  these 
interpolation  functions  do  not  accurately  account  for  the 
physical  behaviour  of  the  flow.  The  errors  generated  by 
these  interpolation  functions  thus  affect  the  solution  in 
the  normal  direction  by  diffusion  and  in  the  downstream 
direction  by  convection.  As  a  result,  the  downstream 
calculation,  particularly  in  the  turbulent  near  wake 
region,  suffers  from  these  errors. 

According  to  the  theory  of  Alber  [62],  the  mean  x- 
directional  velocity  component  has  the  logarithmic 
behaviour  of  the  turbulent  boundary  layer  in  the  turbulent 
inner  near  wake  region  although  it  breaks  down  as  the  flow 
moves  downstream.  In  the  present  study,  the  first  x- 
directional  computational  grid  in  the  wake  region  is 
placed  at  the  starting  point  of  turbulent  inner  wake 
region  (  x"*"  ^  250  )  and  the  first  three  y-directional  u- 
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velocity  calculation  nodes  are  placed  at  the  logarithmic 
layer  (  y"^  ^  60,  150,  250  )  .  Thus,  the  computational 
nodes  used  in  the  interpolation  function  for  the  initial 
portion  of  turbulent  inner  near  wake  region  are  placed  in 
the  logarithmic  layer.  The  discrepancies  between  the 
physical  behaviour  and  interpolation  functions  result  in  a 
considerable,  error,  especially  in  the  evaluation  of  u- 
velocxty.  If  the  square  interpolation  function  is  used  in 
the  evaluation  of  u-velocity,  it  gives  a  higher  u- 
velocity  values  for  the  second  node  and  the  wake 
centerline.  Consequently,  this  higher  u-velocity 
influences  the  generation  term  of  kinetic  energy  in  the 
third  node.  Thus,  the  turbulent  kinetic  energy  is 
underpredicted.  These  errors  can  be  reduced  if  many  finer 
grids  are  placed  near  the  wake  centerline.  However,  there 
is  a  limit  because  the  wall  function  used  in  the  boundary 
layer  calculation  can  only  be  applied  to  the  region  60  < 
y"*"  <  400.  On  the  other  hand,  only  the  laminar  sublayer, 
y"*"  <  10,  from  the  boundary  layer  flow  is  destroyed  in  the 
starting  portion  of  turbulent  inner  wake  region.  In  other 
words,  the  wake  symmetry  condition  satisfies  only  about  a 
laminar  sublayer  thickness  above  the  wake  centerline. 

Thus,  an  accurate  resolution  of  the  symmetry  condition 
along  the  wake  centerline  by  a  simple  interpolation 
function  is  impossible,  and  more  importantly,  the  solution 
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will  be  quite  grid  dependent.  It  is  quite  apparent  that  a 
direct  imposition  of  the  symmetry  condition  (J)(i,2)  =<j>(i,l) 
over  a  distance  of  Ay"*"  in  the  order  of  60  should  cause 
severe  error. 

One  possible  way  of  avoiding  these  errors  will  be  a 
direct  calculation  along  the  wake  centerline  with  the 
inlet  condition  provided  properly.  If  we  consider  the 
calculation  of  asymmetric  wake,  the  whole  domain  should  be 
calculated.  Rhie  and  Chow  [16]  performed  the  whole  domain 
calculation  even  in  the  symmetric  wake  calculation  to 
avoid  numerical  error  from  interpolation  along  the  line  of 
symmetry.  However,  details  of  the  inlet  conditions  of 
transport  quantities  for  the  wake  calculation  immediately 
downstream  of  the  plate  were  not  explained  in  their  paper. 

Since  the  k  -  e  turbulence  model  is  used  in  this 
study,  the  control  volume  which  contains  the  laminar  wake 
region  should  be  excluded  from  the  wake  calculation. 

Since  the  laminar  sublayer  and  buffer  layer  and  a  part  of 
the  logarithmic  layer  in  the  boundary  layer  on  the  body 
are  not  calculated  in  the  wall  function  method,  the 
detailed  destruction  of  these  layers  in  the  wake  can  not 
be  calculated.  Thus  a  universal  function  which  can 
resolve  the  large  streamwise  velocity  gradient  caused  by 
the  destruction  of  these  layers  should  be  used  in  the  wake 


function  region.  The  idea  of  creating  the  wake  function 
is  the  same  as  that  the  wall  function  is  introduced  in  the 
calculation  of  boundary  layer  on  the  body. 
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IV. 4. 2.  Proposed  Wake  Functions. 


In  order  to  adopt  a  wake  function,  we  note  that  Alber 
[62]  obtained  the  mean  velocity  profile  in  the  turbulent 
inner  near  wake  region  using  the  boundary  layer 
approximation  and  similarity  transformation  under  the 
assumption  of  a  linear  distribution  of  eddy  viscosity  in 
the  normal  direction.  The  wake  centerline  velocity  and 
velocity  components  near  the  wake  centerline  in  the 
turbulent  inner  near  wake  region  are 


u 


TO 


where 


-  [  In  g(x+)  -  Y  ]  +  B 

K 

^  [  In  y+  +  Ei(C)  ]  +  B 
K  [  1  -  exp(-C)  3  /  In  g(x+) 


(IV-10) 

(IV-11) 

(IV-12) 


u^o  y'  ^  “  y+/g(x+) 

oo 

Ei(C)  =  /  exp(-t)/t  dt 

and 


g(a)  [  In  g (a)  -  1  ]  =  k2  a 


(IV-13) 
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and  is  the  trailing  edge  friction  velocity  and  y 
(  =  0.5772157  )  is  the  Euler  constant.  The  constants  K,  B 
in  these  equations  were  taken  as  K  =  0.42,  B  “  5.5  for  the 
consistency  with  the  wall  function.  Eq. (IV-12)  was 
obtained  from  Eq. (IV-11)  and  continuity  equation. 

One  may  alternatively  use  the  centerline  velocity 
formula  given  in  Eq. (IV-10)  with  the  correlation  obtained 
from  experiment  by  Andreopoulos  and  Bradshaw  [66]; 
u 

— ^  =  2.02  In  x+  +  0.7  (IV-14) 

to 

Alber's  solution,  Eq. (IV-10),  is  used  as  the  inlet 
condition  for  the  downstream  wake  centerline  u-velocity 
calculation  in  this  study. 

For  the  proper  introduction  of  wake  functions  for  the 
turbulent  quantities  in  the  wake  function  region,  we  may 
consider  Chevray  and  Kovasnay's  data  [63].  They  show  that 
the  turbulent  structure  in  the  logarithmic  layer  from  the 
trailing  edge  of  the  flat  plate  to  the  starting  point  of 
turbulent  inner  near  wake  region  remains  the  same  although 
the  structure  of  laminar  sublayer  and  buffer  layer  evolve 
rapidly  by  means  of  mixing  in  the  laminar  wake  region. 
Therefore,  we  may  use  the  second  y-directional  nodal 
values  of  the  k  and  e  at  the  trailing  edge  of  plate  as  the 
second  y-directional  nodal  values  of  k  and  e  in  the  first 
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control  volume  in  the  wake,  as  shown  in  Fig.IV-2,  without 
loss  of  accuracy.  Obviously,  a  direct  calculation  of 
these  quantities  is  impossible  since  the  calculation  point 
for  kw  and  ew  is  located  immediately  downstream  of  the  the 
trailing  edge  where  the  centerline  u-velocity  can  not  be 
calculated.  We  thus  adopt 


-1/2  2 

k  =  c,,  u„ 

w  |X  to 


e  =  u,^  /  Ky 
w  to 


(IV-15) 

(IV-16) 


as  the  second  nodal  values  of  the  turbulent  kinetic 
energy  and  the  rate  of  its  dissipation  in  the  first 
control  volume  in  the  wake.  Here,  u^^  is  the  trailing  edge 
friction  velocity. 


There  is  no  detailed  experimental  observation  made 
for  the  variation  of  k  and  e  along  the  wake  centerline 
within  our  knowledge.  This  lack  of  inlet  conditions  for  k 
and  e  along  the  wake  centerline  lead  us  to  use  the 
interpolation  function  (  (|>  =  Ay2  +  B  )  in  the  subsequent 
calculation  of  k,  e  along  the  wake  centerline.  However, 
this  does  not  create  error  in  the  calculation  of  k,  e  in 
the  second  node  in  the  wake  because  the  value  of  k  and  e  in 
the  second  node  in  the  wake  are  more  dependent  on  the 
generation  term  of  the  turbulent  kinetic  energy  than  the 
k,  e  values  along  the  wake  centerline.  In  other  words. 
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accurate  resolution  of  the  wake  centerline  velocity  is 
more  important  for  the  accurate  prediction  of  second  nodal 
turbulent  kinetic  energy  than  the  specification  of  k  and  e 
along  the  wake  centerline. 

IV. 4. 3  Enforcing  the  Conservation  of  Mass 

An  important  problem  in  the  calculation  of  wake  by 
the  wall  and  wake  function  method  is  the  calculation  of 
the  second  nodal  v-velocity  component  in  the  wake  region. 
The  difficulty  arises  from  the  fact  that  the  detailed 
destruction  (  mixing  )  of  laminar  sublayer  and  buffer 
layer  and  logarithmic  layer  in  the  first  y-directional 
control  volume  in  the  wake  can  not  be  properly  described 
in  the  wall  function  method  due  to  the  lack  of 
computational  nodes  near  the  wake  centerline.  Obviously, 
this  is  more  severe  in  the  calculation  of  the  intial 
portion  of  turbulent  inner  near  wake.  To  avoid  these 
numerical  difficulties,  we  adopt  PSL  (  Parabolic  Sublayer 
)  like  treatment  [73] .  In  this  study,  the  second  nodal  v- 
velocity  in  the  wake  which  is  denoted  as  V2  in  Fig.IV-2  is 
calculated  by  the  mass  conservation  of  the  first  control 
volume  using  the  following  simple  integration  formula. 

Aye  Aye 

Axc  ^2  ”  J  Uj^  dy  -  J  U2  dy  (IV-17) 

0  0 


where 


^yc  1  13  1 

j  u  dy  =  (  j  u(i,l)  +  u(i,2)  +  ^  u(i,3)  )  Aye 

0 

(IV-18) 


The  accuracy  of  the  integration  formula  is  very 
significant  in  the  wake  calculation  since'  a  small  error  in 
mass  conservation  may  lead  to  a  large  error  in  the 
prediction  of  the  pressure  field  in  the  weak  interaction 
like  present  problem.  It  is  not  easy  to  find  a  good 
integration  formula  which  can  accurately  resolve  the  mass 
flux  of  laminar  sublayer,  buffer  layer  and  a  part  of 
logarithmic  layer  at  the  trailing  edge  of  body  (  x  =  0  ) 
so  that  v^  in  Fig.IV-2  can  be  calculated  by  Eq. (IV-17) . 
Thus  the  Eq. (IV-12)  is  used  only  for  the  specification  of 
the  v-velocity  component  in  the  wake  function  region  (  v^ 
in  Fig.IV-2  )  for  simplicity  and  accuracy.  The  pressure 
in  the  wake  is  calculated  from  the  third  node  in  the  y- 
direction  using  these  v-velocity  components  as  velocity 
boundary  condition.  The  u-velocity  component  in  the 
second  node  is  calculated  using  the  pressure  in  the  third 
node.  The  starting  point  of  calculation  in  the  y- 
direction  in  the  wake  region  is  shown  in  Fig.IV-2. 
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IV. 5.  Results  and  Discussions 

Calculations  were  performed  using  the  partially 
parabolic  form  of  finite  analytic  method  which  is  derived 
in  chapter  II-5-1  and  SIMPLER  numerical  algorithm  outlined 
in  chapter  III-4.  The  derivation  of  pressure  and  pressure 
correction  equation  is  given  in  chapter  III-2.  The  time 
marching  technique  is  employed  for  the  iteration  process 
and  satisfactory  convergence  were  obtained  after  100  time 
marching  with  time  step  At  =  1. 

The  calculation  is  made  for  Reynolds  number  2.48  x 
10®  based  on  the  plate  length  and  free  stream  velocity. 

The  predicted  results  are  presented  and  compared  with 
available  experimental  data  of  Ramaprian  et  al.[64]  who 
measured  the  mean  velocity  and  turbulent  quantities  in  the 
turbulent  inner  near  wake  region. 

The  converged  pressure  distribution  near  the  trailing 
edge  of  the  flat  plate  along  y"*"  ~  150  over  the  plate  and 
wake  centerline  is  shown  in  the  Fig.IV-3.  Fig.IV-3  shows 
the  pressure  near  the  trailing  edge,  either  upstream  or 
downstream,  is  lower  than  the  free  stream  pressure  which 
is  set  equal  to  zero.  The  minimum  pressure  is  located  at 
the  trailing  edge  where  the  fluid  experiences  the  maximum 
acceleration  and  deceleration.  According  to  Alber  [62], 
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the  pressure  gradient  effect  is  mainly  confined  to  a  small 
laminar  wake  region  close  to  the  trailing  edge  where  x"*"  < 
100.  Obviously,  the  present  wall  and  wake  function  method 
together  with  the  k  -  e  turbulence  model  can  not  predict 
the  detailed  interaction  in  this  laminar  wake  region. 
However,  the  sharp  reduction  followed  by  a  rapid  increase 
of  pressure  near  the  trailing  edge  and  its  slow  recovery 
in  the  wake  is  predicted  by  the  present  method  as  shown  in 
Fig.IV-3.  The  prediction  of  pressure  distribution  by 
Patel  and  Chen  [70]  by  two-layer  model  is  given  in  Fig. IV- 
4 .  There  exists  a  small  difference  between  the  present 
prediction  and  the  calculation  by  Patel  and  Chen  [70]  in 
the  region  of  plate.  However,  both  predictions  give 
nearly  the  same  magnitude  of  trailing  edge  interaction. 

The  calculated  result  of  the  skin  friction 
coefficient  is  shown  in  Fig.IV-5  and  compared  with  the 
conventional  flat -plate  correlation.  In  general,  Fig.IV-5 
shows  a  good  agreement  between  the,  predicted  value  and 
simple  flat  plate  correlation.  Since  the  pressure  on  the 
plate  drops  rapidly  near  the  trailing  edge  causing  the 
boundary  layer  to  accelerate,  the  skin  friction 
coefficient  near  the  trailing  edge  is  expected  to 
increase.  This  increase,  even  though  it  is  small,  is  well 
predicted  in  the  present  calculation. 
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Fig.IV-6  shows  the  comparison  of  the  predicted  wake 
centerline  velocity  distribution  with  experimental  data  by 
Ramaprian  et  al.[64]  The  first  computational  node  along 
the  centerline  is  located  approximately  at  x"*"  ~  250.  An 
excellent  agreement  is  obtained  by  the  application  of  wall 
and  wake  function  method.  As  compared  with  the  predicted 
wake  centerline  velocity  distribution  by  Patel  and  Chen 
[70],  which  is  shown  in  Fig.IV-7,  the  present  wake 
function  method  accurately  predicts  the  wake  development 
without  the  detailed  calculation  of  the  laminar  sublayer 
on  the  plate  and  the  laminar  wake. 

Figs . IV-8, IV-9,  and  IV-10  show  the  detail  of  the 
predicted  u-velocity  component,  Reynolds  shear  stress  and 
turbulent  kinetic  energy.  The  comparison  of  the  predicted 
results  with  experimental  data  is  excellent.  As  mentioned 
before,  the  direct  imposition  of  symmetry  condition  or  the 
use  of  interpolation  function  as  the  symmetry  condition 
for  u-velocity  along  the  centerline  may  lead  to 
predictions  of  a  higher  u-velocity  and  lower  turbulent 
kinetic  energy  due  to  the  inadequate  grid  concentration 
near  the  wake  centerline  when  the  wall  function  method  is 
adopted  in  the  calculation  of  boundary  layer  on  the  body. 
On  the  other  hand,  the  present  method  of  combining  the 
wall  and  wake  function  clearly  removes  these  difficulties. 
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It  should  be  mentioned  that  while  the  selection  of 
appropriate  wake  functions  is  important,  an  accurate 
evaluation  of  the  second  nodal  v-velocity  component  is 
also  important  in  achieving  the  accurate  prediction  of  the 
near  wake  turbulent  structure.  We  note  that  the  turbulent 
kinetic  energy  and  turbulent  shear  stress  measurements  of 
Ramaprian  et  al .  [64]  give  slightly  higher  values  than 
those  of  Pot  [65]  in  the  initial  portion  of  turbulent 
inner  near  wake  region  which  can  be  seen  at  the  locations 
X  =  25.4  mm  and  x  =  127  mm. 

A  proper  prediction  of  the  v-velocity  component  in 
the  turbulent  near  wake  region  is  very  siginificant  for 
the  accurate  prediction  of  the  flow  structure  since  the 
amount  of  mixing  in  this  region  is  more  or  less  related  to 
the  v-velocity  component .  The  second  nodal  v-velocity  in 
the  wake  region  calculated  from  the  present  PSL  (parabolic 
sublayer)  like  treatment  can  not  be  confirmed  due  to  the 
lack  of  experimental  data.  A  comparison  with  Alber  theory 
[62]  is  made.  Two  calculations  are  performed  with  a 

+ 

different  grid  distribution  near  the  wake  centerline.  Ay 

c 

+ 

~  120  (  grid  1  ) ,  Ay  ~  200  (  grid  2  )  where  Ay  is  the 

c  c 

size  of  the  first  y-directional  control  volume  in  the  wake 
as  shown  in  Fig.IV-2.  Fig.IV-11  shows  the  comparison  of 
the  present  prediction  with  Alber  theory  [62]  in  the 
turbulent  inner  near  wake  region.  Good  agreement  is 
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obtained  but  the  integral  formula  gives  a  little  grid 
dependency . 

Fig.IV-12  and  Fig.IV-13  indicate  the  prediction  of 
velocity  defect  and  Reynolds  shear  stress  in  the  far  wake 
region  where  the  velocity  defect  w  (  =  1  -  u  )  becomes 
small  compared  with  u-velocity  component.  The  assumption 
of  similarity  of  velocity  and  shear  stress  profiles  leads 
the  half  power  law  for  the  decay  of  the  centerline 
velocity  defect  wq  (  =  1  -  uc  )  and  the  growth  of  the  half 
width  ^  where  w/wq  =*  1/2  )  .  According  to 

asymptotic  theory  [71], 

b  ~  xl/2  ^  wq  ..  x”l/2  (IV-19) 

A  simple  velocity  defect  and  shear  stress 
distribution  can  be  obtained  [69]  when  one  introduces  a 
constant  eddy  viscosity  across  the  wake  in  the  normal 
direction  for  the  momentum  equation. 

w/wo  =  exp(  -  4  ln2  T]2  )  (IV-20) 

_  2 

-uv/w^  =  4  (  7C  ln2  )  1/2  (Vt/u0)  T]  exp  (  -  4  ln2  112  ) 

(IV-21) 

where  T]  =  y/b,  and  Vt/u0  =  0.032  which  has  been 
confirmed  by  Rodi  [74]  from  a  survey  of  several  sets  of 
experimental  data.  Fig.IV-12  and  Fig.IV-13  show  the 
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comparison  between  present  prediction  and  asymptotic 
theory,  Eq. (IV-20)  and  Eq. (IV-21) .  The  velocity  defect, 
in  general,  is  predicted  well  except  near  the  edge  of  the 
wake  where  the  computed  profile  underpredicted  the 
asymptotic  value.  This  may  be  due  to  the  fact  that  the 
turbulence  model  can  not  properly  simulate  the  high 
intermittance  of  turbulent  and  laminar  flow  occurring  at 
the  edge  of  the  wake .  On  the  other  hand,  the  turbulent 
shear  stress  is  underpredicted.  This  is  a  well  known 
defect  of  the  current  k  -  e  turbulence  model  which  was 
already  observed  in  the  previous  calculations  [67],  [69], 

and  [70]  .  An  improvement  in  the  turbulence  modelling  is 
required  in  order  to  achieve  a  better  agreement  in  the 
prediction  of  turbulent  structure  in  the  far  wake  region 
with  the  experiment.  For  example,  Chen  and  Singh  [75] 
showed  a  better  prediction  can  be  achieved  with  the  k  -  E 
model  based  on  the  two  turbulence  scale  concept  where  the 
first  scale  is  based  on  the  turbulent  kinetic  energy,  k, 
and  its  dissipation  rate,  E,  to  characterize  the  large, 
energy-containing  eddies  while  the  second  scale,  the 
Kolmogorov  scale,  is  based  on  the  dissipation  reate,  E,  and 
the  kinematic  viscosity,  V,  to  characterize  the  small, 
energy-dissipating  eddies. 

A  more  sophisticated  analysis  of  wake  functions  and 
related  calculation  method  is  needed  for  the  application 


of  the  present  method  to  the  calculation  of  wakes  behind 
bodies  of  streamline  curvature  for  which  the  pressure 
gradient  effect  are  important. 
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lYxi^CAnclusion 

A  new  calculation  method  is  presented  in  the  present 
study  for  an  accurate  prediction  of  turbulent  wake  flow 
behind  a  flat  plate.  The  method  proposes  a  wake  function 
that  models  the  laminar  and  turbuent  wake  immediately 
downstream  of  the  trailing  edge  of  the  plate  from  x"*"  =  0 
to  approximately  x'*’  =  250.  The  combination  of  the  wake 
function  for  the  near  wake  region  and  the  wall  function 
for  the  near  wall  region  on  the  plate  provides  a  means  for 
accurate  prediction  of  turbulent  wake  flow  and  eliminates 
excess  computational  grids  required  near  the  wall  and  the 
trailing  edge  and  the  uncertainty  of  the  turbulence  model. 
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CHAPTER  V 

LAMINAR  AND  TURBULENT  FLOW  PAST  A  FINITE  FLAT  PLATE 

V.l  Introduction 

The  classical  problem  of  incompressible  laminar  and 
turbulent  flow  past  a  finite  flat  plate  has  been  studied 
by  numerous  investigators.  It  is  a  problem  that  not  only 
provides  much  fundamental  understanding  of  basic  feature 
of  leading  and  trailing  edge  interaction,  boundary  layer 
on  the  plate  and  wake  development,  but  also  of 
considerable  importance  in  diverse  engineering 
applications  such  as  thin  airfoil  theory. 

When  the  Reynolds  number  based  on  the  plate  length 
and  the  free  stream  velocity  is  large,  the  laminar 
boundary  layer  on  the  plate  is  described  by  the  well-known 
Blasius  solution.  However,  the  Blasius  solution,  which  is 
based  on  the  boundary  layer  theory  of  Prandtl,  does  not 
properly  describe  the  flow  field  near  the  leading  and 
trailing  edges  of  plate  due  to  the  breakdown  of 
assumptions  made  in  the  boundary  layer  theory. 

Due  to  its  geometric  simplicity,  the  laminar  flow 
near  the  trailing  edge  of  plate  has  been  subject  of 
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various  numerical  and  theoretical  investigations.  The 
study  of  this  problem  provides  us  the  basic  understanding 
of  inviscid-viscous  interaction  at  the  trailing  edge  and 
near  wake  development .  According  to  the  triple  deck 
theory  of  Stewartson  [76]  or  Messiter  [77],  for  infinite 
Reynolds  number  flow,  the  laminar  boundary  layer  on  the 
plate  develops  into  Goldstein's  near-wake  through  a  small 
transition  region  around  the  trailing  edge  which  is  known 
as  the  triple-deck  region.  This  triple  deck  region  which 
is  embedded  inside  the  boundary  layer  near  the  trailing 
edge  arises  in  order  to  avoid  the  ‘ .  jularity  of  the 
boundary  layer  equations  at  the  trailing  edge.  Melnik  and 
Chow  [78]  and  Veldman  and  van  de  Vooren  [79]  have  obtained 
the  numerical  solutions  of  the  triple  deck  equations . 

The  breakdown  of  boundary  layer  theory  for  describing 
the  flow  field  near  the  trailing  edge  of  plate  leads  many 
investigators  to  use  alternative  approaches.  Two 
approaches  are  commonly  used:  the  inviscid-viscous 
interaction  method  [80] ;  and  the  numerical  solutions  of 
fully  elliptic  or  partially  parabolic  form  of  equations 
[81,82].  Recently,  Chen  and  Patel  [81]  have  obtained  the 
numerical  solutions  of  this  problem  using  the  elliptic 
form  of  finite  analytic  numerical  scheme.  The  grid 
dependence  and  convergence  tests,  the  influence  of  the 
size  of  solution  domain  on  the  solution  were  investigated 
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in  detail.  The  predicted  results  were  compared  with 
previous  results  by  interacting  boundary  layer  theory  and 
by  partially  parabolic  numerical  solution.  One  valuable 
observation  made  in  this  study  is  that  the  solution  domain 
should  be  larger  than  that  used  in  the  previous  studies  in 
order  to  capture  the  whole  inviscid-viscous  interaction 
and  far  wake  development  and  to  obtain  domain  independent 
solution.  The  same  observation  is  also  reported  in  Rhie 
[83]  in  the  calculation  of  laminar  flow  past  a  thin 
airfoil . 

However,  most  of  the  previous  studies  are  confined  to 
the  analysis  of  the  flow  field  near  the  trailing  edge. 

Thus  the  leading  edge  interaction  and  the  initial 
development  of  boundary  layer  near  the  leading  edge  were 
not  considered.  In  order  to  remove  the  dependency  of  the 
solution  on  the  location  of  inlet  solution  boundary  and  to 
obtain  whole  flow  field  past  a  plate,  the  calculations 
should  be  started  at  the  upstream  of  the  plate  with  inlet 
conditions  provided  by  uniform  flow  conditions.  The 
earlier  numerical  solutions  of  the  full  Navier-Stokes 
equations  obtained  by  Dennis  and  Dunwoody  [84]  are  limited 
to  relatively  low  Reynolds  number  flow.  Xu  and  Chen  [82] 
has  recently  obtained  a  numerical  solution  of  complete 
laminar  flow  field  past  a  finite  flat  plate  using  a  novel 
nine-point  finite  analytic  numerical  scheme.  The 
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predicted  results  of  pressure  distribution  on  the  plate 
and  along  the  wake  centerline,  the  wake  centerline 
velocity  development  and  the  friction  coefficient 
distribution  near  the  trailing  edge  were  presented  with 
comparision  with  previous  studies. 

In  the  present  study,  the  classical  problem  of 
laminar  and  turbulent  flow  past  a  finite  plate  is 
reconsidered.  The  laminar  flow  past  a  finite  flat  plate 
is  first  solved  using  both  partially  parabolic  form  of 
finite  analytic  numerical  method  given  in  chapter  II-5-1 
and  elliptic  form  of  finite  volume  method  with  modified 
TEACH  code.  This  practice  will  not  only  give  a  critial 
evaluation  of  suitability  of  the  partially  parabolic  form 
of  Navier-Stokes  equations  in  the  prediction  of  the  strong 
leading  edge  interaction  but  also  provide  us  the 
comparisions  of  results  by  both  methods  which  has  not  been 
made  before.  Then  the  numerical  solutions  of  turbulent 
flow  past  a  finite  flat  plate  is  obtained  by  solving  the 
Reynolds  averaged  Navier-Stokes  equations  with  the  k  -  e 
turbulence  model  and  the  wall  function  method.  The  finite 
volume  method  is  used  in  this  calculation.  This  practice 
will  show  how  the  introduction  of  turbulence  model 
influences  the  overall  solution.  The  predicted  results 
are  compared  with  available  experimental  data  of  Ramaprian 
et  al.  [64]  . 
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V.2.  Numerical  Grids  and  Solution  Domains 

The  numerical  grids  are  generated  using  the  method 
outlined  in  chapter  II-4-1  and  the  partial  view  of 
generated  grids  are  shown  in  Fig.II-3  for  both  laminar  and 
turbulent  calculations.  Table-3  shows  the  details  of 
solution  domain  and  number  of  computational  grids  used  in 
the  present  calculations .  The  solution  domain  in  normal 
direction  for  laminar  flow  calculation  is  chosen  to  be 
~  12L  in  which  Chen  and  Patel  [81]  obtained  the  domain 
independent  solution.  The  solution  domain  in  the  other 
directions  is  also  made  large  enough  to  capture  the  whole 
leading  and  trailing  edge  interactions.  The  generated 
numerical  grids  are  treated  as  control  volume  lines  and 
the  computational  grid  lines  are  placed  at  the  center  of 
control  volume  lines. 

_ Boundary  Conditions 

The  upstream  boundary  conditions  are  specified  by 
following  uniform  flow  conditions; 

U  -  1,  =  0,  k  =  k^^,  E  =  (V-I) 

where  k.  and  e,  is  specified  as  follows 
in  in 

k^^  =  1.5  (  Tu 


(V-2) 
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''in^  '  ^In 


(V-3) 


and  the  turbulent  intensity,  Tu  =  0.5%  and  length 
scale,  Ij^j^  =  0.03L  are  used  in  the  present  calculations. 


For  laminar  flow  calculations,  the  no-slip  condition 
is  specified  on  the  plate.  The  wall  function  approach  of 
Launder  and  Spalding  [43]  is  used  for  turbulent  flow 
calculation.  The  other  boundary  conditions  are  same  as 
those  reported  in  chapter  II-3-2. 


V.4  Numerical  Method 

The  finite  analytic  computer  code  used  in  the 
previous  chapter  is  modified  for  the  laminar  flow 
calculation.  The  well-known  TEACH  code  [86]  is  used  for 
finite  volume  calculations  with  subtantial  modification. 
The  solution  procedures  outlined  in  chapter  III-4  are 
implemented  to  TEACH  code  which  includes  the  change  of 
numerical  algorithm,  from  SIMPLE  to  SIMPLER  and  the  usage 
of  parabolic  marching  technique  with  global  pressure 
calculation  procedure  instead  of  fully  elliptic 
calculation  procedure.  Details  of  TEACH  code  are  well- 
documented  in  reference  [86]  and  are  not  explained  here. 
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V.5  Results  and  Dlscuaalons 

V.5.1  Laminar  Flow  Past  A  Finite  Flat  Plate 

Calculations  are  performed  for  Re  ”  10^  using  both 
partially  parabolic  finite  analytic  method  outlined  in 
chapter  II-5-1  and  elliptic  finite  volume  method  given  in 
reference  [86] .  The  same  numerical  grids  are  employed  for 
both  calculations.  Of  135  x-directional  numerical  grids, 
35  grids  are  placed  upstream  of  plate,  60  grids  are  placed 
on  the  plate  and  40  grids  are  placed  in  the  wake  region. 
The  first  y-directional  u-calculation  point  over  the  plate 

_3 

is  made  Ay  =  0.35x10  and  the  minimum  grid  size  at  the 

leading  edge  is  made  (Ax)j^  =«  0.897x10  and  the  minimum 

_2 

grid  size  at  the  trailing  edge  is  made  (Ax)^  =  0.2x10 
Time  marching  technique  is  employed  for  finite  analytic 
calculation  with  time  step  At  =  1 .  The  relaxation  factors 
used  in  the  finite  volume  calculation  are  =  0.7 

and  ttp  =  0.5.  Satisfactory  convergence  was  obtained  after 
500  iterations  for  the  calculation  by  finite  volume  method 

-3 

(  RES  =  0.6617x10  )  and  1000  iterations  for  the 

calculation  by  partially  parabolic  finite  analytic  method 

-3 

(  RES  =  0.2891x10  )  where  RES  denotes  sum  of  mass 

residuals  in  pressure  correction  equation  divided  by  inlet 
flow  rate.  It  was  observed  that  the  velocity  components 
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and  pressure  field  are  settled  down  within  200  iterations 
and  the  rest  computational  efforts  are  devoted  to 
reduction  of  the  mass  residuals  to  a  certain  degree. 

The  predicted  pressure  distribution  on  the  plate  and 
along  the  wake  centerline  is  shown  in  Fig.V-1.  As  shown 
in  the  figure,  the  magnitude  of  leading  edge  interaction 
is  lager  than  the  magnitude  of  trailing  edge  interaction. 
However,  the  interacting  region  near  the  leading  edge  is 
smaller  than  the  interacting  region  near  the  trailing 
edge.  The  comparision  of  predicted  pressure  distributions 
by  two  different  numerical  schemes  shows  that  both  methods 
predicted  nearly  same  pressure  distribution  except  near 
the  leading  edge.  The  supurious  pressure  drop  near  the 
leading  edge  predicted  by  the  present  partially  parabolic 
numerical  scheme  is  not  observed  in  the  present  elliptic 
calculation  or  in  the  nine-point  elliptic  calculation  by 
Xu  and  Chen  [82] .  This  fact  shows  that  the  partially 
parabolic  assumption  may  not  be  appropriate  for  the 
simulation  of  strong  leading  edge  interaction.  However, 
it  is  also  shown  that  the  partially  parabolic  assumption 
is  valid  for  the  prediction  of  most  of  the  boundary  layer 
on  the  plate,  the  trailing  edge  interaction  and  wake 
development . 
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The  enlarged  view  of  pressure  distribution  near  the 
trailing  edge  is  shown  in  Fig.V-2.  The  magnitude  of 
trailing  iteraction  predicted  by  finite  volume  method  (  ~ 
-0.01563  )  is  slightly  higher  than  that  predicted  by 
partially  parabolic  finite  analytic  method  (  ~  -0.01474  ). 
The  origin  of  this  difference  is  not  clearly  understood. 
For  reference,  the  trailing  edge  pressure  predicted  by 
Chen  and  Patel  (81 ]  using  elliptic  finite  analytic  method 
is  reported  as  -0.0152. 

Fig.V-3  shows  the  distribution  of  friction 
coefficient  near  the  trailing  edge  by  both  methods.  The 
difference  in  friction  coefficient  near  the  trailing  edge 
by  both  methods  is  consistent  with  the  predicted  pressure 
distribution  in  this  region.  The  stronger  interaction 
cause  stronger  acceleration  near  the  trailing  edge  and 
larger  penetration  to  the  upstream  of  the  plate.  This 
figure  also  clearly  shows  that  the  Blasius  solution  which 
is  based  on  the  boundary  layer  theory  is  not  adequate  for 
describing  the  flow  field  near  the  trailing  edge. 

The  development  of  wake  centerline  velocity  is  also 
affected  by  the  trailing  edge  interaction  when  the 
accelerated  boundary  layer  at  the  trailing  edge  is 
destroyed  in  the  wake  with  adverse  pressure  gradient . 
Fig.V-4  shows  the  predicted  wake  centerline  velocity 
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profiles  in  the  wake  region  by  two  different  numerical 
methods.  Little  difference  is  evident,  which  is 
consistent  with  the  pressure  distribution  in  this  region. 
There  may  exist  some  effects  which  are  caused  by  using 
different  numerical  treatments  near  the  wake  centerline. 

In  the  finite  volume  method,  the  wake  centerline  velocity 
do  not  affect  the  second  nodal  velocity  in  the  wake  by 
placing  the  control  volume  surface  along  the  wake 
centerline  and  applying  the  symmetry  condition  and  zero 
normal  velocity  condition.  In  finite  analytic  method,  the 
wake  centerline  velocity  affects  the  second  nodal 
velocity.  These  practice  may  result  in  smaller  wake 
centerline  velocity  prediction  by  finite  analytic  method. 
However,  this  effect  seems  to  be  very  small  since  the  wake 
centerline  velocity  prediction  by  the  present  elliptic 
finite  volume  method  more  agree  well  with  those  reported 
in  Chen  and  Patel  [81]  who  used  elliptic  finite  analytic 
method.  It  is  also  noted  that  these  two  elliptic 
calculations  have  nearly  same  magnitude  of  trailing  edge 
interaction.  Although  there  may  exist  some  numerical 
errors,  it  may  be  concluded  that  the  fully  elliptic  form 
of  discretization  equation  should  be  used  in  order  to 
accurately  capture  the  leading  and  trailing  edge 
interactions  and  thereby  accurately  predict  the  wake 
development . 
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However,  the  differences  of  predictions  by  two 
different  methods  are  very  small  and  are  confined  to  a 
relatively  small  region.  These  effects  are  not 
significant  in  the  practical  calculations  when  the 
interaction  is  not  weak  as  in  the  present  problem. 

V.5.2  Turbulent  Flow  Past  A  Finite  Flat  Plate 

g 

Calculations  are  performed  for  Re  =  2.48x10  , 

corresponding  to  experiments  of  Ramaprian  et  al .  [64], 

using  finite  volume  method  with  the  standard  k  -  £ 

turbulence  model  and  the  conventional  wall  function 

method.  Of  120  x-directional  numerical  grids,  30  grids 

are  placed  upstream  of  plate,  55  grids  are  placed  on  the 

plate  and  next  35  grids  are  placed  in  the  wake  region. 

The  first  y-directional  u-calculation  point  over  the  plate 

is  made  Ay  =  0.164x10  which  correspondins  to  y  -  15  at 

the  trailing  edge.  The  minimum  grid  size  at  the  leading 

—4 

edge  is  made  (Ax)  =  0.323x10  and  the  minimum  grid  size 

at  the  trailing  edge  is  made  (^x) ^  =  0.4x10  .  The 

relaxation  factors  used  in  this  calculation  are  a  =  a  * 

u  V 

=  Og  =  0.8  and  Op  =  0.5.  Satisfactory  convergence  was 
obtained  after  500  iterations.  (  RES  =  0.4203x10”^  ) 

Fig.V-5  and  Fig.V-6  show  the  predicted  pressure 
distribution  on  the  plate  and  along  the  wake  centerline. 
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There  exists  a  sharp  peak  at  the  leading  edge  and  a  sharp 
drop  at  the  trailing  edge  followed  by  slower  recovery  in 
the  wake  region.  The  predicted  leading  and  trailing  edge 
interactions  in  turbulent  flow  are  much  weaker  than  those 
in  laminar  flow  by  the  introduction  of  turbulence  model 
and  insufficient  grid  refinement  in  these  regions.  The 
grid  refinement  in  these  regions  is  restricted  since  the  k 
-  e  turbulence  model  can  only  be  applied  to  fully  turbulent 
region.  Thus  the  detailed  interactions  near  the  leading 
edge  and  in  the  laminar  wake  region  near  the  trailing  edge 
can  not  be  captured.  Patel  and  Chen  [70]  explained  that 
the  slower  r<*c'  very  of  wake  centerline  pressure  to  ambient 
value  is  t‘  lated  to  the  isotropic  eddy  viscosity 
as sump cion  in  the  k  -  E  turbulence  model; 


9p  ^  2^  ^ 

dx  3  dx 


(v-4) 


The  predicted  friction  velocity  distribution  on  the 
plate  is  shown  in  Fig.V-7  with  a  comparision  of  the  simple 
flat  plate  correlation  given  in  Reference  [87] .  The 
strange  transition  near  the  leading  edge  region  is  due  to 
the  improper  modelling  of  transition  by  k  -  e  turbulence 
model  which  is  also  observed  in  Rhie  [83]  in  the 
calculation  of  turbulent  flow  past  a  thin  airfoil.  There 
exist  a  substantial  difference  between  the  prediction  and 
correlation  even  in  fully  turbulent  region.  It  was  found 
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that  the  underprediction  of  friction  velocity  originates 
from  use  of  the  wall  function  method.  The  first  u- 
velocity  calculation  point  over  the  plate  is  placed  around 
y^  ~  15  at  the  trailing  edge  to  provide  proper  calculation 
in  the  wake  region.  Thus  all  the  first  calculation  points 
on  the  plate  are  placed  in  the  buffer  layer.  However,  the 
wall  function  region  employed  in  the  actual  calculation  is 
divided  into  only  two  regions,  laminar  flow  region  (  y^  < 
11.63  )  and  fully  turbulent  region  (  y^  >  11.63  ).  Thus 
the  wall  shear  stress  on  the  plate  is  calculated  by  way  of 
logarithmic  law  even  though  the  first  u-calculation  point 
is  placed  in  the  buffer  layer.  As  a  result,  the  friction 
velocity  is  underpredicted  even  though  the  velocity  is 
overpredicted.  The  friction  velocity  in  the  turbulent 
region  is  corrected  using  the  third  nodal  u-velocity 
component  which  lie  in  logarithmic  layer  and  the  law  of 
the  wall.  The  corrected  skin  friction  velocity  is 
consistent  with  the  overprediction  of  velocity  at  the 
trailing  edge  which  is  shown  in  Fig.V-9.  The  amount  of 
overprediction  in  friction  velocity  is  not  so  serious  as 
shown  in  the  figure.  It  may  be  suggested  that  the 
predictions  will  be  improved  if  the  wall  shear  stress  in 
the  wall  function  method  is  evaluated  using  the  third 
nodal  velocity  component  in  the  practical  calculation.  It 
was  also  found  that  the  prediction  of  friction  velocity 
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was  slightly  improved  when  the  first  grid  over  the  plate 
was  slightly  increased  to  y^  ~  30.  However,  this  practice 
was  abandoned  due  to  the  very  poor  prediction  in  the  wake 
region  which  was  explained  in  detail  in  the  previous 
chapter.  It  was  also  observed  that  the  prediction  was 
extremely  sensitive  to  the  numerical  grid  distribution 
near  the  wall  region  which  was  not  observed  in  the 
calculation  of  previous  chapter  in  which  two  point  wall 
function  method  was  adopted.  This  sensitivity  of 
prediction  to  the  location  of  the  first  grid  point  is  due 
to  the  Couette-flow  type  wall  function  method  usually  used 
in  the  finite  volume  calculations.  It  is  noted  that  the 
two-point  wall  function  method  can  not  be  used  in  the 
laminar  flow  region  near  the  leading  edge.  The  poor 
prediction  of  friction  velocity  can  be  improved  with  more 
advanced  turbulence  model  which  can  properly  describe  the 
flow  field  near  the  viscosity-affected  near  wall  region. 

Figs.V-8,  V-9  and  V-10  show  the  predicted  wake 
centerline  velocity,  mean  velocity  ditribution  and 
turbulent  kinetic  energy  distribution  in  the  near  wake 
region.  The  predicted  results  are  only  fairly  good  due  to 
the  improper  calculation  of  the  boundary  layer  on  the 
plate.  The  overprediction  of  velocity  profile  at  the 
trailing  edge  results  in  higher  velocity  profiles  and 
smaller  kinetic  energy  distribution  in  the  near  wake 
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region.  However,  these  predictions  are  not  much  deviated 
from  the  experimental  data  as  shown  in  these  figures.  It 
is  shown  that  as  the  flow  approaches  downstream,  the  flow 
becomes  to  lose  the  memory  of  boundary  layer  on  the  plate 
and  the  predictions  agree  better  with  the  experimental 
data . 


Although  only  fairly  good  predictions  were  achived 
due  to  improper  modelling  of  transition  and  near  wall 
turbulence,  the  present  study  achieves  the  prediction  of 
whole  turbulent  flow  past  a  finite  flat  plate.  The 
predictions  can  be  improved  if  a  suitable  turbulence  model 
is  used.  However,  there  does  not  exist  a  general 
turbulence  model  which  can  properly  describe  the 
transition  of  laminar  to  turbulent  and  the  flow  field  in 
the  viscosity-affected  near  wall  region. 
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CHAPTER  VI 

TURBULENT  FLOW  PAST  AXISYMMETRIC  BODIES 

VI. 1  Introduction 

In  order  to  design  an  optimal  ship  body  and  ship 
propelling  system,  it  is  very  important  to  have  a 
fundamental  understanding  and  to  be  able  to  make  an 
accurate  prediction  of  fluid  flow  past  various  ship 
bodies .  The  accurate  prediction  of  turbulent  flow  near 
the  ship  stern  region  is  particularly  important  since  many 
propellers  and  appendages  are  located  inside  the  ship 
stern  boundary  layers.  The  flow  evolution  in  this  ship 
stern  region  is  characterized  by  a  rapid  growth  of 
boundary  layer,  a  strong  viscous-inviscid  interaction  and 
a  general  reduction  in  the  level  of  turbulence  by 
streamline  curvature  and  pressure  gradient  [88]  which  is 
quite  different  from  that  in  the  hull  region.  These  flow 
features  show  that  the  thin  boundary  layer  approximation 
is  not  adequate  for  describing  the  flow  field  around  the 
ship  stern  region  although  it  may  be  suitable  for 
describing  the  flow  field  around  the  ship  hull  region. 
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A  review  of  the  state-of-the-art  of  available 
experimental  and  numerical  investigations  of  turbulent 
flow  past  axisymmetric  bodies  was  recently  made  by  Patel 
and  Chen  [89] .  According  to  this  detailed  survey,  many 
reliable  experimental  informations  are  available  at 
present  which  can  be  used  to  verify  the  accuracy  of 
numerical  calculations.  Chevry  [90]  made  measurements  of 
mean  flow  and  turbulent  qpaantities  in  the  wake  of  6  to  1 
spheroid.  These  data  provide  quite  detailed  informations 
on  the  near  wake  evolved  from  the  a  separated  flow  just 
ahead  of  the  tail  region  of  spheroid.  The  experiment  of 
Patel  et  al.  [88]  was  conducted  on  the  same  model  but  the 
separation  was  eliminated  by  means  of  a  short  conical  tail 
attachment.  They  measured  the  development  of  thick 
boundary  layer  on  this  modified  spheroid.  Patel  and  Lee 
[91]  made  extensive  measurements  on  the  boundary  layer  and 
near  wake  development  of  turbulent  flow  past  a  low-drag 
body  (  F-57  body  ) .  Since  above  two  measurements  are  made 
normal  to  the  body  surface  using  boundary  layer 
coordinates,  it  is  difficult  to  use  these  data  to  compare 
with  practical  calculation  if  numerical  grids  are  not 
generated  in  the  same  manner.  The  extensive  experimental 
data  measured  by  Huang  et  al.  [92,93,94]  for  four 
different  afterbodies,  afterbody-1,  afterbody-2, 
afterbody-3,  afterbody-5,  provide  the  detailed 
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informations  on  the  thick  boundary  layer  and  the  near  wake 
development .  These  four  bodies  have  the  same  streamlined 
forebody  and  parallel  middle  body,  but  have  different 
afterbodies.  All  of  these  bodies  have  an  inflection  point 
and  quite  dramatic  changes  in  surface  curvature  which 
induce  a  strong  favorable  and  adverse  pressure  gradient 
over  the  stern  region  and  propellor  hub  region.  In  the 
case  of  afterbody-3,  experimental  data  indicate  that  there 
is  a  small  separation  bubble  around  the  inflection  point. 
These  data  have  been  widely  used  in  the  previous  studies 
[89,95,96,97,98,99]  for  the  evaluation  of  calculation 
methods  and  turbulence  models.  These  data  are  also  used 
in  the  present  study  to  compare  with  the  calculated 
results . 

Many  numerical  efforts  have  been  made  in  the  past 
decade  to  compute  this  complex  flows  near  the  tail  region 
of  the  axisymmetric  bodies.  These  efforts  have  involved 
with  different  approximations  made  in  the  governing 
equations,  with  different  turbulence  models  and  with 
different  calculation  methods.  According  to  the  survey  by 
Patel  and  Chen  [89],  the  previous  calculation  methods  can 
be  broadly  classified  into  three  categories,  the  non¬ 
interactive  numerical  solution  of  thin  layer  or  boundary 
layer  equation  with  the  pressure  field  given  by 
experimental  data  [91],  the  viscous-inviscid  interaction 
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method  [100,101,102]  and  the  global  numerical  solutions  of 
partially  parabolic  or  elliptic  ec[uations  [89,98],  Among 
these  methods,  the  global  numerical  solution  method  is  the 
most  attractive  since  this  method  does  not  require  the 
experimental  informations  of  pressure  field  or  another 
numerical  solution  of  inviscid  flow.  It  is  interesting  to 
note  that  the  recently  developed  calculation  methods 
outlined  in  chapter-I  have  not  been  applied  to  this 
problem  except  the  works  of  Patel  and  Chen  [89]  .  The 
calculation  methods  used  by  Muraoka  [97,103]  or  Markatos 
[98]  are  not  general  and  can  only  be  applicable  when  the 
numerical  grids  are  generated  in  a  special  algebraic 
manner.  It  is  also  noted  that  most  of  previous 
calculations  start  at  the  middle  of  body  to  avoid  the 
numerical  difficulties  associated  with  the  calculation 
near  the  leading  edge . 

In  the  present  study,  numerical  solutions  of 
partially  parabolic  or  fully  elliptic  form  of  Reynolds- 
averaged  Navier-Stokes  equations  in  numerically  generated, 
nonorthogonal  coordinates  are  obtained  through  the 
calculation  procedure  outlined  in  chapter-II  and  chapter- 
ill.  The  complex,  physical  geometry  is  resolved  by  the 
use  of  nonorthogonal,  body-fitted  coordinates  which  are 
generated  through  the  solutions  of  Poisson  equations  with 
appropriate  grid  control  functions.  The  governing 


118 


equations  are  written  in  the  transformed  coordinates  by 
the  partial  transformation  where  the  original  orthogonal 
velocity  components  are  left  as  dependent  variables. 

These  transformed  governing  equations  are  discretized 
using  the  finite  analytic  method  by  Chen  and  Chen  [3] 
which  was  already  explained  in  chapter-II.  The  finite 
volume  method  is  also  applied  to  some  problems  in  order  to 
provide  comparisions  of  two  different  numerical  schemes  in 
the  prediction  of  turbulent  boundary  layer  and  wake  of 
axisymmetric  bodies .  The  one-velocity  staggered  grid 
technique  is  employed  for  all  the  calculations  presented 
in  this  chapter.  The  derivation  of  pressure  and  pressure 
correction  equations  in  this  grid  configuration  is  given 
in  chapter  III-3-1. 

Since  there  are  strong  pressure  gradient  and 
streamline  curvature  effects  on  the  flow  near  the  stern 
region,  some  considerations  should  be  made  for  the 
selection  of  til^rbulence  model.  The  zero  equation  models 
of  Cebeci-Smith  model  or  Baldwin-Lomax  model  and  the  one 
equation  model  have  been  used  in  the  earlier  calculations 
[91,95,96,101]  with  some  modifications.  However,  the  k  -  e 
turbulence  model  with  wall  function  method  has  been  the 
most  widely  employed  in  the  recent  calculations 
[89,97,102,103].  The  defect  of  the  k  -  e  turbulence  model 
in  the  calculation  of  flow  involving  pressure  gradient 
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effect  is  well  explained  by  Rodi  and  Scheuerer  [54] .  The 
poor  performance  of  the  current  k  -  e  turbulence  model  is 
caused  by  inaccurate  modelling  of  the  e  -  equation.  For 
example,  the  length  scale  determined  using  the  e  -  equation 
rises  steeper  near  the  wall  region  in  the  adverse  pressure 
gradient  region  while  the  experimental  data  shows  that  the 
length  scale  is  virtually  independent  of  the  presure 
gradient.  By  this  reason,  the  turbulence  models  using  the 
empirical  length  scale  specification  near  the  wall  region 
yield  much  better  predictions  for  the  adverse  pressure 
gradient  boundary  layers  than  does  the  current  k  -  e 
turbulence  model.  Rodi  and  Scheuerer  [54]  suggest  use  of 
Hanjaliac  and  Launder's  model  [104] .  On  the  other  hand, 
Chen  and  Patel  [56]  developed  a  two-layer  turbulence  model 
to  consider  the  pressure  gradient  and  streamline  curvature 
effects  and  flow  separation.  Their  model  has  been 
successfully  applied  to  the  prediction  of  turbulent  flow 
past  axisymmetric  bodies . 

In  the  present  study,  some  calculations  are  performed 
using  the  conventional  k  -  e  turbulence  model  and  wall 
function  method  in  which  the  pressure  gradient  effect  is 
embodied  in  the  wall  function  while  some  calculations  are 
performed  using  the  two-layer  turbulence  model  of  Chen  and 
Patel  [56] . 
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VI  .,2 _ Two-Layer  Model 

In  the  two-layer  model,  the  calculation  region  is 
divided  into  two  regions;  a  near  wall  region  and  an  outer 
region,  as  shown  in  the  Fig. VI-1.  The  boundary  between 
the  near  wall  region  and  the  outer  region  is  placed  around 
y+  ~  100.  The  near  wall  region  includes  the  laminar 
sublayer,  buffer  layer  and  a  part  of  logarithmic  layer 
ajacient  to  the  body.  In  this  region,  calculations  are 
carried  all  the  way  to  the  wall  using  the  one  equation 
turbulence  model  where  the  rate  of  turbulent  kinetic 
energy  dissipation  e  and  the  turbulent  eddy  viscosity  is 
specified  by  algebraic  relations  to  consider  the  wall 
proximity  viscous  effects.  The  outer  region  is  the  upper 
region  of  the  near  wall  region  over  the  body  and  the  whole 
wake  region.  The  standard  k  -  e  turbulence  model  is  used 
to  simulate  this  fully  turbulent  region. 

Three  one  equation  turbulence  models  such  as  Hassid 
and  Poreh's  model  [48],  Norris  and  Reynolds's  model  [49], 
Wolfshtein's  model  [47]  were  tested  in  the  course  of 
present  investigation.  However,  the  Wolfshtein  model  is 
chosen  in  the  present  study  by  the  simple  reason  that  the 
other  two  models  slightly  overpredict  the  wall  shear 
stress  compared  with  the  Wolfshtein  model. 
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In  the  Wolfshtein  model,  the  turbulent  eddy  viscosity 

and  the  rate  of  turbulent  kinetic  energy  dissipation  e 

is  specified  by  a  simple  algebraic  form  as  following; 

,  3/2 


e  = 


,  1/2  , 

Vt  =  k  1^ 


(VI-1) 

(VI-2) 


where  the  length  scales  1^,  1^  contain  the  viscosity 
damping  effect  in  the  near  wall  region  in  terms  of 
turbulent  Reynolds  number  Rj^. 


1|^  =  Cl  n  [  1  -  exp(  ’  Rjc  /  )  1 

Ig  *  Cl  n  [  1  -  exp(  -  /  Ag  )  ] 


(VI-3) 

(VI-4) 


where 


Rj^  =  Re  n  k^^^. 


0.09,  cj^  =  K  c 


-3/4 


and  n  is  the  normal  distance  from  the  wall. 


These  constants  are  evaluated  by  Chen  and  Patel  [56] 
and  somewhat  different  from  those  reported  in  Wolfshtein 
[47]  or  Yap  [58]  . 


The  numerical  grids  are  generated  using  the  grid 
generation  technique  outlined  in  chapter  II-4-2.  The 
generated  numerical  grids  for  four  different  afterbodies 
are  shown  in  Fig. I 1-4.  The  detailed  informations  on  the 
solution  domains  are  given  in  Table  4  and  Table  5.  These 
solution  domains  are  large  enough  to  capture  the  whole 
viscous-inviscid  interactions  and  wake  development. 


In  the  half  body  calculations  using  the  wall  function 
method,  the  upstream  inlet  conditions  were  specified  by 
following  simple  correlations. 


~  In  y^  +  B 

K 

-1/2  2  y 

%  '"t  ^  ^  ■  5  ^ 

3/4  3/2  ,  ^ 

c^  k  /  Ky 


(VI-5) 

(VI-6) 

(VI-7) 


where  y^  =  Re  u  y,  K  is  von-Karman  constant,  K  =  0.42 

and  B  =  5.5.  The  boundary  layer  thickness  5  and  friction 

velocity  u  =  (T  /p)^/2  are  provided  by  experimental  data. 
X  w 


The  wall  function  method  employed  in  the  present 
study  is  a  little  different  from  the  conventional  wall 


function  method  usually  used  in  the  control  volume  method. 
In  the  present  study,  it  is  assumed  that  at  least  two  u- 
calculation  nodes  are  located  in  the  logarithmic  region. 


The  friction  velocity  u^  is  first  determined  by  iterative 
method  using  the  velocity  components  in  the  third  node 
which  is  obtained  from  the  previous  iteration  and 
following  the  generalized  form  of  the  law  of  the  wall  [4] 
in  which  the  pressure  gradient  effects  on  the  flow  in  the 
wall  region  are  taken  into  consideration. 


(  1  +  n+  ) 1/2  -  1 

=  -  {  In  [  - ; - 1 — 772 -  ] 

K  A^  (  1  +  n+  )l/2  +  1 

+2  [  (  1  +  A_  n+  )l/2  -1]  }+B+3.7A  (VI-8) 

X  p 


where  n"^  =  Re  u_n  is  dimensionless,  normal  distance 
X 

3 

from  the  surface,  A  =  -  Ap  /  {  Re  u  )  is  the 

P  X 

dimensionless  pressure  gradient  on  the  body  surface,  is 
the  dimensionless  shear  stress  gradient  and  is  assumed  to 

be  7  A  q  is  the  magnitude  of  velocity  or  q  =  (  u  +  v^ 
1/2 

)  ,  K  =  0.418  is  the  von-Karman  constant,  and  B  =  5.45. 


Then  the  magnitude  of  velocity  in  the  second  node  q2 
can  be  calculated  using  this  friction  velocity  and  Eq. (VI- 
8)  while  the  velocity  components  parallel  to  ^  -  grid  line 
are  calculated  by  the  following  relations . 
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Vill 


r  (  b,  +  b. 


q|  =  q^ 

^3 

where  the  geometric  coefficients  J  and  b 

2  2 

in  Eqs. (11-58)  and  (11-59)  and  ~  ^  ’ 


j 


(VI-9) 

(VI-10) 


are  given 


The  u,  V,  k,  e  values  at  the  second  node  are 
evaluated  using  the  friction  velocity  and  following 
relations . 


^2  =  -r=^  q; 

V^ii  ^ 

(VI-11) 

''2  ~  ~r^  % 

\9ll  ^ 

(VI-12) 

-1/2  2 

>^2  =  %  ^ 

(VI-13) 
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62  =  /  Ky 

(VI-14) 

The  details  of  this  two-point  wall  function  method  is 
given  in  Chen  and  Patel  [4] .  The  advantage  of  this  two- 
point  wall  function  method  is  that  the  sensitivity  of  the 
solution  to  the  location  of  the  first  grid  point,  which  is 
found  in  other  Couette-flow  type  wall  function  method 
usually  used  in  the  control  volume  method,  is  removed. 


In  the  two-layer  calculations,  the  upstream  u- 
velocity  boundary  condition  near  the  wall  is  specified  by 
following  Spalding's  velocity  law; 
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The  inlet  boundary  conditions  for  turbulent 
quantities  in  the  inner  layer  are  specified  by  following 
correlations  given  in  Patel  and  Chen  [70]  while  turbulent 
quantities  in  the  outer  layer  is  specified  by  Eq. (IV-6) 
and  Eq. (IV-7) . 


+  +2 
k  =  0.05  (  y  ) 

=  1.25  +  0.325  (  y"^  -  5  ) 

=  4.5  +  (  y'*’  -  15  )  /  37.5 

=  3.3 


y'*’  <  5  (VI-16) 

5  <  y"^  <  15 
15  <  y"*”  <  60 
60  <  y"*"  <  120 


and 

e'^  =  0.1  y"^  /  120 
=  1  /  Ky"^ 


y'*’  <  12  (VII-17) 

y""  >  12 


The  dimensionless  quantities  in  these  equations  are 

,  +  ^  2  +  4  + 

defined  as  k  =  k  /  u^  ,  e  =  e  /  (  Re  u^  )  ,  y  =  Re  u^  y 

and  u^  =  u  /  u„ . 

X 


The  other  boundary  conditions  are  same  as  those 


reported  in  chapter  II-3-2. 


VI. 5.1  Prediction  of  Turbulent  Flow  Past  Axisymmetric 
Bodies  using  Wall  Function  Method 

Calculations  are  performed  for  three  different 

axisymmetric  bodies  of  DTNSRDC,  afterbody-1,  afterbody-2 

and  afterbody-5.  Both  finite  analytic  method  and  finite 

volume  method  are  used  to  solve  Reynolds-averaged  Navier- 

Stokes  equations  with  the  k  -  £  turbulence  model  and  two- 

point  wall  function  method.  The  same  numerical  grids  and 

boundary  conditions  are  employed  for  both  calculations. 

The  time  marching  technique  is  employed  for  finite 

analytic  calculations  with  time  step  At  =  1  and  pressure 

relaxation  factor  a  *  0.2.  The  relaxation  factors  used  in 

P 

the  finite  volume  calculations  are  a  =■  a  =  0.7,  a,  = 

u  V  '  k  e 

=  0.8  and  a  =0.5. 

P 

Figure  VI-2  shows  the  convergence  history  of  both 
calculations  for  the  afterbody-1.  The  calculation  by 
finite  volume  method  converges  faster  than  that  by  finite 
analytic  method  in  the  earlier  stage.  However,  the 
calculation  by  finite  volume  method  has  a  limit  in 
convergence.  This  limited  convergence  originates  from  use 
of  one  velocity  staggered  grid  method  which  is  also 
reported  in  Shyy  et  al.  [24].  As  will  be  shown  in  the 
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next  chapter,  use  of  two  velocities  staggered  grid  method 
eliminates  this  deficiency.  The  calculation  by  finite 
analytic  method  monotonically  converges  without  a  limit 
although  the  same  one  velocity  staggered  method  is  used. 

-3 

Satisfactory  convergence,  RES  <  2x10  ,  was  obtained 

within  100  iterations  for  both  calculations  where  RES 
denotes  the  sum  of  mass  residuals  in  the  pressure 
correction  equation  divided  by  the  inlet  flow  rate. 
However,  the  pressure  and  velocity  fields  are  settled  down 
within  50  iterations  as  shown  in  Figs. VI-3  to  VI-8.  These 
figures  show  that  convergences  of  boundary  layer  and  wake 
development  by  finite  analytic  method  are  a  little  faster 
than  those  by  finite  volume  method.  This  is  due  to  the 
fact  that  the  relaxation  factors  in  the  finite  volume 
method  is  same  in  the  whole  computational  domain  while 
relaxation  factors  in  the  finite  analytic  method  are 
changed  according  to  the  size  of  numerical  grid  and  the 
magnitude  of  convecting  velocity  components  which  is 
another  novel  nature  of  the  finite  analytic  method. 

Figs. VI-9  to  VI-11  show  the  converged  pressure 

2 

distributions,  Cp  =  2  (p-p  )  /  (pu  ),  on  the  body  surface 

oo 

and  along  the  wake  centerline  for  afterbody-1,  afterbody- 
2,  and  afterbody-5.  As  shown  in  the  figures,  the 
agreement  between  the  predictions  and  measured  data  is 
excellent.  The  differences  in  the  prediction  by  two 
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different  methods  are  small  and  are  confined  to  only  the 
region  of  inlet,  propellor  hub  and  trailing  edge.  The 
dramatic  change  of  pressure  in  the  stern  and  propeller  hub 
region  is  well  predicted.  The  pressure  is  recovered  in 
the  wake  at  a  distance  roughly  0.3L  from  the  trailing  edge 
of  the  body.  Since  the  solution  domain  of  the  present 
study  extends  far  beyond  this  distance,  the  whole  viscous- 
inviscid  interaction  is  completely  captured. 

Figs. VI-12  to  VI-14  show  the  comparision  between 
predicted  and  measured  friction  velocity.  Fairly  good 
agreement  between  predictions  and  measurments  was  obtained 
in  the  stern  region.  The  slight  disagreement  in  the 
upstream  region  may  due  to  the  simple  imposition  of  inlet 
conditions.  The  results  can  be  improved  if  the  accurate 
inlet  conditions  of  velocity  and  turbulent  quantities  are 
provided  by  experimental  data.  The  finite  volume  method 
consistently  overpredicts  the  friction  velocity  in  all 
calcultions.  As  shown  in  the  figures,  there  is  steep  drop 
of  friction  velocity  in  the  adverse  pressure  gradient 
region  and  the  flow  is  very  near  to  separate  at  x  =  0.93L 
for  afterbody-5,  x  =  0.97L  for  afterbody-2  but  no  flow 
separation  is  predicted. 

F:ugs. VI-15  to  VI-17  show  the  radial  pressure 


distribution  in  the  stern  and  near  wake  region. 
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Satisfactory  agreements  between  the  prediction  and 
measurements  are  made  although  they  are  not  as  good  as 
those  of  pressure  distribution  on  the  body  surface.  The 
small  discrepancies  may  be  partly  due  to  the 
dif ficultities  in  the  measurements.  The  differences  in 
prediction  between  two  methods  are  small  and  are  again 
confined  to  the  regions  of  propellor  hub  and  trailing 
edge.  We  note  that  the  viscous-inviscid  interaction 
extends  radially  to  0.35L  from  the  surface  of  the  body 
near  the  trailing  edge.  This  fact  shows  the  failure  of 
thin  boundary  layer  approximation  for  the  present  problem 
and  some  careful  considerations  should  be  made  for  the 
selection  of  radial  direction  solution  domain.  In  the 
present  study,  the  computational  domain  extends  radially 
to  y^  ~  1 . IL  as  shown  in  Table  4 . 

Comparisions  of  the  predictions  with  the  measured 
data  for  both  axial  and  radial  velocity  components  is  seen 
in  Figs. VI-18  to  VI-20.  The  agreements  between  measured 
data  and  the  predictions  by  finite  analytic  method  is 
excellent.  The  finite  volume  method  slightly  overpredicts 
the  axial  velocity  component,  especially  near  the  tail 
region  and  in  the  near  wake  region.  One  may  notice  the 
thickening  of  boundary  layer  in  the  tail  and  in  the  near 
wake  region  which  is  consistent  with  the  pressure 
variation  in  these  regions. 
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-he  turbulent  kinetic  energy  in  the  tail  region  and 
in  the  near  wake  region  is  overpredicted  as  shown  in  the 
Figs. VI-21  to  VI-23.  The  amount  of  overprediction  is  more 
severe  in  the  prediction  of  afterbody-5  and  afterbody-2 
than  that  in  the  afterbody-1.  The  differences  in 
prediction  by  two  methods  are  small  and  consistent  with 
the  predicted  velocity  profiles.  It  is  noted  that  the 
predictions  of  turbulent  kinetic  energy  made  by  Huang  and 
Chang  [99]  agrees  quite  well  with  experimental  data.  They 
used  a  different  near  wall  treatment  instead  of  the  wall 
function  method.  As  will  be  shown  in  the  next  section, 
the  use  of  two-layer  model  slightly  improve  the 
prediction.  From  these  observations,  the  overprediction 
of  turbulent  kinetic  energy  in  the  tail  and  near  wake 
regions  seems  to  be  originated  from  the  inadequacy  of  k  -  e 
turbulence  model  in  the  prediction  of  fluid  flow  involving 
pressure  gradient  and  streamline  curvature  effects.  The 
other  predictions  made  by  Patel  and  Chen  [89] ,  Muraoka 
[97,103]  using  the  k  -  e  turbulence  model  and  the  wall 
function  method  show  the  similiar  trend  as  in  the  present 
prediction . 

Overall  agreements  between  measurements  and 
predictions  made  using  both  methods  are  encouraging.  For 
most  of  the  stern  region  and  near  wake  region,  the 
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pressure  and  velocity  components  are  accurately  predicted. 
But  some  considerations  must  be  made  on  the  selection  of 
turbulence  model  especially  on  the  near  wall  treatment  in 
order  to  accuately  predict  the  turbulent  quantities  for 
the  flow  involving  pressure  gradient  and  streamline 
curvature  effects.  Although  predictions  by  finite 
analytic  method  are  in  favor  for  all  the  calculations,  the 
differences  in  the  predictions  by  two  different  numerical 
schemes  are  small.  It  shows  that  the  selection  of 
numerical  scheme  is  not  important  for  calculating  of 
boundary  layer  type  flow  where  the  numerical  false 
diffusion  effects  are  very  small.  The  overall  success  of 
present  calculations  are  mainly  due  to  the  use  of  two- 
point  wall  function  method.  However,  the  use  of  two-point 
wall  function  method  is  limited  and  can  not  be  confidently 
applied  to  the  flow  involving  separation.  Calculations 
are  also  made  using  the  partially  parabolic  form  of  finite 
analytic  method  although  it  is  not  presented  here.  The 
differences  in  predictions  between  partially  parabolic 
calculations  and  elliptic  calculations  were  found  to  be 
negligible.  This  fact  shows  that  partially  parabolic 
assumption  is  a  good  approximation  if  there  is  no  flow 
separation . 
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VI. 5. 2  Prediction  of  Turbulent  Flow  Past  Axisynunetric 
Bodies  by  Two-Layer  Model 

Calculations  are  performed  for  four  different 
axisynunetric  bodies  of  DTNSRDC,  afterbody-1,  afterbody-2, 
afterbody-3  and  afterbody-5.  The  partially  parabolic  form 
of  finite  analytic  method  explained  in  chapter  II-5-3  is 
employed  to  solve  Reynolds-averaged  Navier-Stokes 
equations  with  two-layer  model.  The  FLARE  approximation 
where  the  ^-directional  convecting  velocity  (  in  Eq. (II- 
84)  )  is  zero  if  it  is  less  than  zero  is  employed  in  order 
to  properly  handle  the  flow  separation  ocurring  during  the 
iteration  process.  The  introduction  of  FLARE  scheme  is 
come  from  the  observation  made  by  Ramakrishnan  and  Rubin 
[105]  .  They  found  that  the  FLARE  scheme  was  the  most 
stable  for  the  calculation  of  small  separated,  high 
Reynolds  number  flow.  The  time  marching  technique  is 
employed  for  iteration  process  with  time  step  At  =  0.002. 
Satisfactory  convergence  wc.  >  obtained  after  1000  time 
marching . 

Figs.IV-24  to  IV-27  show  the  converged  pressure 

2 

distributions,  Cp  =  2  (p-p  )  /  (pU  ),  on  the  body  surface  and 
along  the  wake  centerline.  As  shown  in  the  figures,  the 
agreements  between  the  measured  data  and  the  predicted 
results  are  fairly  good.  The  pressure  is  a  little 


sensitive  to  the  shape  of  the  body  surface  as  compared 
with  predictions  using  the  wall  function  method.  It  is 
due  to  the  fact  that  pressure  is  calculated  all  the  way  to 
the  wall  in  the  two-layer  calculations  while  pressure  on 
the  body  surface  is  extrapolated  from  the  wall  function 
region  to  the  body  surface  using  a  simple  linear  function 
in  the  calculations  using  the  wall  function  method.  It  is 
noted  that  the  predictions  for  afterbody-3  are  possible  by 
use  of  the  two-layer  model. 

Figs.IV-28  to  VI-31  show  the  comparisions  between  the 
predicted  and  measured  friction  velocity.  Athough  fairly 
good  agreements  are  made  in  the  hull  region,  systematic 
overpredictions  are  made  in  the  tail  region.  The 
overprediction  of  friction  velocity  in  the  tail  region  may 
be  partly  due  to  the  numerical  difficulties  in  the 
generation  of  fine  numerical  grid  near  the  wall  in  the 
tail  region  of  the  body  and  may  be  partly  due  to  the  fact 
that  even  the  one  equation  turbulence  model  tends  to 
overestimate  the  velocity  in  the  near  wall  region 
indicating  a  reduction  of  turbulence  length  scale  in  the 
tail  region.  It  is  noted  that  flow  separation  was  not 
predicted  in  the  present  calculation  for  afterbody-3  in 
which  experimental  data  indicate  there  exist  a  small 
separation  bubble  around  the  inflection  point.  Figs. VI-32 
to  VI-33  show  the  predicted  friction  velocity  for 
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afterbody-3  and  afterbody-5  by  Chen  and  Patel  [56]  using 
the  same  two-layer  model  and  the  nine-point  finite 
analytic  method.  There  exist  substantial  differences 
between  two  predictions.  The  origin  of  these  differences 
can  not  be  clearly  explained  since  the  details  of 
pressure-velocity  coupling  technique  are  not  reported  in 
Chen  and  Patel  [56] .  However,  the  present  predictions 
does  not  show  the  strange  overshoots  which  is  observed  in 
the  predictions  by  Chen  and  Patel  [56] . 

Fairly  good  agreements  between  the  predictions  and 
the  measured  data  were  obtained  for  both  the  axial  and 
radial  velocity  components  which  can  be  seen  in  Figs. VI-34 
to  VI-37  except  the  axial  velocity  component  is  slightly 
overpredicted  in  the  tail  region  as  explained  before. 

The  overprediction  of  turbulent  kinetic  energy  in  the 
tail  region  and  in  the  near  wake  region  in  the 
calculations  using  the  wall  function  method  is  slightly 
improved  by  using  the  two-layer  model  as  shown  in  the 
Figs. VI-38  to  VI-40.  This  observation  shows  that  the 
defect  of  the  k  -  e  turbulence  model  in  the  calculation  of 
turbulent  boundary  layer  and  wake  involving  pressure 
gradient  and  curvature  effects  is  a  little  relieved  by  the 
introduction  of  one  equation  turbulence  model  for  the 
calculation  of  flow  near  the  wall  region. 
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Overall  agreements  between  measurements  and 
predictions  made  by  the  present  finite  analytic  method 
with  both  wall  function  method  and  two-layer  model  are 
satisfactory.  For  most  of  the  stern  region  and  the  near 
wake  region,  velocity  components  and  turbulent  quantities 
are  accurately  predicted  except  that  the  axial  velocity 
component  is  overpredicted  in  the  two-layer  calculation  in 
the  tail  region  of  body.  The  origin  of  this 
overprediction  should  be  clearly  investigated.  It  was 
found  that  the  inlet  condition  severely  influence  the 
overall  solution.  Thus  the  proper  ways  of  specifying  the 
inlet  condition  should  be  sought.  When  it  is  considered 
the  fact  that  the  two-layer  calculation  not  only  provides 
the  numerical  solutions  of  laminar  sublayer,  buffer  layer 
and  a  part  of  logarithmic  layer  without  excessive 
numerical  efforts  but  also  the  two-layer  model  can  be 
applied  to  separated  flow  region  to  which  the  wall 
function  method  can  not  be  confidently  applied,  the  use  of 
two-layer  model  is  highly  recommended. 
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CHAPTER  VII 

TURBULENT  FLOW  PAST  FINITE  AXISYMMETRIC  BODIES 

VTT.l  Introduction 

Incompressible  turbulent  flow  past  axisymmetric 
bodies  of  finite  length  has  been  the  subject  of  numerous 
investigations,  as  it  is  of  major  importance  in 
aerodynamics  and  hydrodynamics .  In  most  of  previous 
studies,  the  calculations  start  at  the  middle  of  body  with 
inlet  conditions  provided  by  curve  fits  of  experimental 
data  [97],  simple  flat  plate  correlations  [56,89]  or 
numerical  solutions  of  boundary  layer  equation  up  to  the 
inlet  location  [103]  .  However,  these  calculations  do  not 
provide  the  solution  of  leading  edge  interaction  and 
initial  development  of  boundary  layer  near  the  leading 
edge.  The  resulting  solutions  are  dependent  on  the  inlet 
conditions  as  shown  in  the  previous  chapter  or  the 
location  of  inlet  boundary  [81] . 

In  order  to  avoid  above  deficiencies  in  the  half  body 
calculations  and  to  obtain  the  solution  of  the  whole  flow 
field,  the  calculation  should  start  far  upstream  of  the 
body  with  inlet  conditions  provided  by  unform  flow 
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condition.  However,  numerical  solutions  of  incompressible 
turbulent  flow  past  axisymmetric  bodies,  including  the 
leading  edge  interaction,  boundary  layer  development  on 
the  body,  the  trailing  edge  interaction  and  the  wake 
development,  is  rarely  seen  in  the  literature.  For  two 
dimensional  flow  situations,  Rhie  and  Chow  [16]  obtained 
the  numerical  solutions  of  turbulent  flow  past  airfoils 
with  or  without  angle  of  attack  using  the  k  -  £  turbulence 
model  and  the  wall  function  method.  The  C-type,  body- 
fitted  coordinates  were  introduced  to  handle  the  complex 
physical  geometry  and  the  SIMPLE  numerical  algorithm  was 
employed  with  non-staggered  grid  arrangement.  Detailed 
numerical  solutions  of  laminar  and  turbulent  flow  past  a 
cicular  cylinder  were  obtained  by  Majumdar  and  Rodi  [55] . 
They  used  an  orthogonal  coordinate  system  with  staggered 
grid  arrangement.  The  k  -  e  turbulence  model  was  also 
employed  in  this  study.  An  improvement  of  prediction  by 
using  the  zero  equation  turbulence  model  in  the  near  wall 
region  instead  of  the  wall  function  method  is  reported. 

Since  the  physical  solution  domain  has  a  complex 
geometry  which  requires  the  use  of  nonorthogonal,  body- 
fitted  coordinate  system,  some  considerations  should  be 
made  about  the  choice  of  calculation  methods  in  this 
coordinate  system.  Among  the  several  ca'' culations  methods 
outlined  in  chapter  I,  the  two  velocities  staggered  grid 
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method  developed  by  Maliska  and  Raithby  [30]  keeps  the 
novel  nature  of  staggered  grid  arrangement  without  any 
numerical  problems  encountered  in  the  use  of  one  velocity 
staggered  grid  method  although  the  use  of  this  grid 
configuration  requires  a  little  more  computer  storage  and 
computing  time  by  storing  both  velocity  components  at  each 
control  volume  surface.  It  is  reported  in  Maliska  and 
Raithby  [30]  that  the  computational  efforts  using  this  two 
velocities  staggered  grid  method  is  not  so  grave  as  it 
looks  by  its  relatively  earlier  convergence.  This  grid 
configuration  allows  to  use  H-type  grid  as  well  as  C-type 
grid  to  which  the  one  velocity  staggered  grid  method  may 
not  be  applicable.  Therefore  the  two  velocities  staggered 
grid  method  is  employed  for  all  the  calculation  presented 
in  this  chapter  with  some  modifications  made  from  the 
original  method  by  Maliska  and  Raithby  [30] . 

In  the  initial  phase  of  present  study,  some 
computational  experiments  were  made  using  the  finite 
analytic  method  and  the  two-point  wall  function  method 
outlined  in  the  previous  chapter.  However,  this  practice 
was  abandoned  since  the  two-point  wall  function  method  did 
not  properly  describe  the  laminar  flow  near  the  leading 
edge  even  in  the  simple  flat  plate  calculation.  It  was 
also  found  that  the  general  law  of  the  wall,  Eq.(VI-8), 
had  singularity  in  the  region  of  favorable  pressure 
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gradient  giving  a  negative  argument  in  the  logarithmic 
function.  Although  this  two-point  wall  function  method 
has  a  salient  feature  of  removing  the  sensitivity  of  the 
first  grid  location  on  the  solution,  the  application  of 
this  method  to  general  flow  situation  is  limited.  After 
testing  other  possibilities,  it  was  decided  to  use  the 
same  numerical  scheme  and  the  near  wall  treatment  as  those 
reported  in  Rhie  and  Chow  [16]  or  Majumdar  and  Rodi  [55] . 
These  authors  used  the  finite  volume  method  and  the 
Couette  flow  type  wall  function  method  by  Launder  and 
Spalding  [43]  .  However,  the  numerical  solutions  using 
this  wall  function  method  usually  suffer  the  sensitivity 
of  the  solution  on  the  numerical  grid  and  the 
overprediction  of  the  axial  velocity  component  in  the 
adverse  pressure  gradient  region  due  to  the  neglection  of 
the  pressure  gradient  effect  in  the  wall  function. 

In  the  present  study,  numerical  solutions  of 
incompressible  turbulent  flow  past  axisymmetric  bodies  are 
obtained  using  the  k  -  e  turbulence  model  and  the  wall 
function  method.  The  complex,  physical  solution  domain  is 
resolved  by  the  use  of  numerically  generated, 
nonorthogonal,  body-fitted  coordinates.  The  governing 
equations  are  written  in  the  transformed  domain  with  the 
cylindrical  veolcity  components  as  the  dependent 
variables.  The  transformed  governing  equations  are 
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discretized  using  the  finite  volume  method  with  hybrid 
numerical  scheme  as  outlined  in  chapter  II-6.  As 
mentioned  before,  the  two  velocities  staggered  grid 
method  is  adopted  and  the  derivation  of  pressure  and 
pressure  correction  equations  in  this  grid  configuration 
is  given  in  chapter  III-3-2. 

Calculation  are  performed  for  four  axisymmetric 
bodies  of  DTNSRDC  using  the  calculation  procedure  outlined 
in  chapter  III-4.  The  numerical  results  are  compared  with 
the  experimental  data. 

Yll.2^  numerical  Grids  and  Solution  Domains 

In  the  initial  phase  of  present  study,  use  of  C-type 
numerical  grid  was  seriously  considered  due  to  the 
roundness  of  forebody  shape.  The  introduction  of  two 
velocity  staggered  grid  method  should  be  understood  within 
this  context  since  the  use  of  one  velocity  staggered  grid 
method  in  C-type  grid  may  cause  numerical  instability. 
However,  the  use  of  C-type  grid  was  abandoned  due  to  the 
difficulties  in  the  specification  of  boundary  conditions. 
The  numerical  difficulties  related  with  the  use  of  C-type 
grid  are  well  reported  in  Rhie  [83]  and  Chen  and  Patel 
[106] .  For  example,  Chen  and  Patel  [106]  first  performed 
calculation  using  H-type  grid.  Then  the  boundary 
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conditions  for  calculation  by  C-type  grid  are  specified  by 
these  calculated  results.  In  order  to  avoid  these 
difficulties,  it  was  decided  to  change  the  forebody  shape 
up  to  =  0.3L  and  to  use  H-type  grid.  The  forebody 
shape  was  changed  to  orgive  shape  using  following  simple 
equations; 


R  =  R. 


-  -  (  X  -  )  ) 


(VII-1) 


where  Rjjj^x  maximum  radius  of  body  and  R^^  is 

defined  as  following  relation. 


^i  ^^max 


(VII-2) 


As  will  be  seen  later,  the  change  of  forebody  shape 
does  not  seriously  influence  the  develoment  of  boundary 
layer  in  the  region  of  middle  body  and  tail  since  these 
DTNSRDC  bodies  have  long  parallel  middle  bodies. 


The  numerical  grids  are  generated  using  the  method 
outlined  in  chapter  II-4-2.  The  partial  view  of  generated 
numerical  grids  are  shown  in  Fig.II-5.  The  x-directional 
numerical  grids  on  the  body  are  distributed  uniformly 
while  those  in  the  upstream  of  body  and  in  the  wake  region 
are  stretched  with  expansion  ratio  1.2  ~  1.3.  Table  6 
shows  the  detailed  information  on  the  flow  condition,  grid 
and  solution  domains  used  in  the  present  calculation. 

These  solution  domains  are  large  enough  to  capture  the 
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whole  leading  and  trailing  edge  interactions.  Of  the  151 
x-directional  numerical  nodes,  20  nodes  are  placed 
upstream  of  body,  101  nodes  are  placed  on  the  body  and  the 
rest  30  nodes  are  placed  in  the  wake. 


VII-3 _ Boundary  conditions 

The  inlet  boundary  conditions  are  specified  by 
following  uniform  flow  conditions. 


u  =  1,  V  =  0,  k  =  k.  ,  e  =  e. 

'  X  '  in'  in 


where  k^.  ^  and  is  specified  as  follows 


(VII-3) 


in 


k^^  =  1.5  (  Tu  ) 


'in  -  =n  ''in^  '  ^In 


(VII-4) 

(VII-5) 


and  the  turbulent  intensity,  Tu  =  0.5%  and  length 
scale,  1^^  =  O.OOIL  are  usea  in  the  present  calculations. 


In  the  wall  function  region,  the  wall  function 
approach  of  Launder  and  Spalding  [43]  is  used.  It  is 
based  on  the  one-dimensional  Couette  flow  analysis.  The 
wall  shear  stress,  is  expressed  as  a  fuxiction  of  the 
second  nodal  velocity  component  parallel  to  the  wall, 
as  follows; 


=  -  X,  U 


w 


P 


(VII-6) 
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where  the  coefficient,  X  is  determined  from  the  two- 

w 

part  universal  velocity  profile. 

X  =  /  5n  if  <  11.63  (VII-7) 

w 

.  1/4  1/2  +  + 

=  p  C|j^  kp  K/  ln(  E  yp  )  if  yp  >  11.63  (VII-8) 

In  these  equations,  K  =  0.418  is  the  von-karman 
constant  and  E  =  0.9793  is  the  roughness  parameter,  5n  is 
the  normal  distance  from  the  wall  and  y^  is  defined  as 
follows; 

+  1/4  1/2  c 

Yp  =  P  kp  5n  /  p  (VII-9) 

Since  the  dependent  variables  of  momentum  equations 
are  cylindrical  velocity  components,  the  shear  force,  = 
should  be  expressed  in  terms  of  its  cylindrical 
components.  The  shear  force,  acting  along  Up  direction 
is  decomposed  into  two  components  and  along  x  and 

y  directions  repectively  as  follows; 

T  =  -  X  A  [  (  1  -  nl  nl  )  u  -  nl  n2  V  ] (VII-10) 

wx  w  w  p  P 

T  -  X  A  [  (  1  -  n2  n2  )  V  -  nl  n2  u  ] (VII-11) 

wy  w  w  p  p 

where  nl  and  n2  are  the  components  of  the  unit  vector 
normal  to  the  wall  along  the  x  and  y  directions 
repectively. 
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ni  =  -  (  /  Vy  )  „ 

n2  =  (  /  Vy  )„ 

and  the  geometric  coefficient  Y 


2  2 


(VII-12) 

(VII-13) 


The  equations  for  turbulent  kinetic  energy,  k  and  its 
dissipation  rate,  e  also  need  a  special  treatment  for  the 
near  wall  cell.  In  the  equations  for  k,  the  diffusion 
flux  through  the  wall  surface  is  set  equal  to  zero  while 
the  source  terms  are  modified  as  follows. 


The  generation  and  dissipation  term  of  turbulent 
kinetic  energy  equation  is  approximated  as  follows  using 
the  Couette  flow  assumption; 


r  G  dv  =  |T„I  |U„|  AV  /  5n 
)Av  ”  P 

f’  3/4  3/2  + 

e  dv  *  Xp  Up  Av  /  8n 


where  U  is  given  by 

Sr 


p  ->fp 


if  Yp  <  11 • 63 


if  yp  >  11.63 


(VII-14) 

(VII-15) 

(VII-16) 

(VII-17) 


The  value  of  e  in  the  near  wall  cell  is  obtained  from 
the  local  value  of  k  using  the  concept  of  local 
equilibrium  between  production  and  dissipation  of 
turbulent  kinetic  energy  as  follows; 
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ep  =  /  (  K  5n  )  (VII-18) 

Since  the  details  of  this  wall  function  method  is 
well-documented  in  Launder  and  Spalding  [43]  or  Peric 
[18],  only  a  brief  description  of  this  method  is  given  in 
the  present  study. 

The  other  boundary  conditions  are  same  as  those 
reported  in  chapter  II-3  and  not  explained  here. 


Calculations  are  performed  for  four  axisymmetric 
bodies  using  the  finite  volume  method  outlined  in  chapter 
II-6.  The  relaxation  factors  used  in  the  present 
calculations  are  a  =a  =0.7,  a,  =a=0.8  and  a  = 

V  jc  &  p 

0.5. 


Fig.VII-1  shows  the  convergence  history  of 

calculations  for  afterbody-1.  The  total  mass  residual 

decreases  rapidly  for  the  first  fifty  iterations.  However, 

the  rate  of  convergence  slows  down  after  fifty  iterations. 

It  is  also  noted  that  the  calculation  using  the  coarse 

grid  gives  a  faster  convergence.  Satisfactory 

-4 

convergence,  RES  <  2x10  ,  was  obtained  within  300 

iterations  although  computations  are  continued  to  500 
iterations  where  RES  denotes  sum  of  mass  residuals  in  the 


146 


pressure  correction  equation  divided  by  inlet  flow  rate. 

In  general,  the  velocity  components  and  pressure  field 
were  settled  down  within  100  iterations  and  the  rest 
computational  efforts  are  devoted  to  the  development  of 
far  wake. 

Figs.VIl-2  and  VII-3  show  the  grid  dependence  tests 
for  pressure  and  friction  velocity.  They  show  that 
present  calculations  are  generally  grid  independent  except 
for  a  small  region  near  the  leading  and  trailing  edges. 

It  should  be  noted  that  the  first  three  y-directional 
computational  grids  near  the  wall  were  fixed  for  both 
calculations  to  avoid  errors  caused  by  the  near  wall 
treatment . 

Figs.VII-4  to  VII-7  show  the  converged  pressure 

2 

distributions,  Cp  =  2 (p-p  ) / (pu  ),  on  the  body  surface  and 
along  the  wake  centerline  for  four  different  afterbodies. 
The  leading  and  trailing  edge  interactions  as  well  as  the 
dramatic  change  of  pressure  in  the  regions  of  stern  and 
propellor  hub  are  well  predicted.  The  agreement  between 
the  predictions  and  measured  data  is  also  fairly  good 
except  pressure  is  overpredicted  in  the  tail  region  of 
afterbody-3.  This  overprediction  is  the  result  of  the 
inadequacy  of  the  present  k  -  e  turbulence  model  with  the 
present  wall  function  method  for  describing  the  flow  field 
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in  this  region.  It  is  noted  that  although  the  magnitude 
of  leading  edge  interaction  is  large,  the  interacting  zone 
is  confined  to  a  small  region.  The  pressure  is  recovered 
in  the  wake  in  a  distance  roughly  0.3L  from  the  trailing 
edge  of  the  body.  Since  the  solution  domain  of  the 
present  study  extends  far  beyond  this  distance,  the  whole 
viscous-inviscid  interaction  is  completely  captured. 

Figs.VII-8  to  VII-11  show  the  predicted  friction 
velocity  distribution  with  comparision  with  measured  data. 
There  exist  a  strange  transition  near  the  leading  edge 
which  is  originated  from  inproper  modelling  of  transition 
by  present  k  -  e  turbulence  model.  The  friction  velocity 
is  underpredicted  in  the  region  of  middle  body  in  spite  of 
velocity  components  in  this  region  are  somewhat  well 
predicted  as  will  be  shown  in  the  next  figures.  As  well 
explained  in  chapter  V,  the  origin  of  this  discrepancy  is 
come  from  the  usage  of  the  wall  function  method  in  the 
fact  that  although  the  first  normal  calculation  point  is 
placed  in  the  buffer  layer  (  y^  ~  20  ) ,  the  friction 
velocity  is  calculated  by  the  logarithmic  law.  The 
increase  of  the  size  of  the  first  grid  would  improve  the 
prediction  of  friction  velocity.  However,  this  practice 
was  abandoned  due  to  the  very  poor  overall  predictions 
especially  in  tail  region  of  the  body.  From  these 
observations,  it  is  suggested  that  the  wall  function 
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region  should  be  at  least  divided  into  three  regions, 
laminar  sublayer,  buffer  layer  and  logarithmic  layer  as 
reported  in  Lee  [107] .  The  use  of  more  elaborated  wall 
functions  by  Launder  and  Chieng  [108]  will  also  improve 
the  predictions.  However,  these  practices  are  not 
performed  in  the  present  investigation.  A  subtantial 
overprediction  of  friction  velocity  is  observed  in  the 
tail  region  of  body,  especially  in  the  predictions  for 
afterbody-3  and  afterbody-5.  These  overprediction  is  due 
to  the  neglect ion  of  pressure  gradient  effect  both  in  the 
turbulence  model  and  in  the  wall  functions.  However,  as 
mentioned  before,  the  use  of  general  law  of  the  wall, 
Eq.(Vl-8),  leads  to  singularity  in  the  region  of  favorable 
pressure  gradient.  It  is  noted  that  flow  separation  is 
not  predicted  for  afterbody-3  although  experimental  data 
i  idicate  that  there  exist  a  small  separation  bubble  around 
the  inflection  point.  It  is  believed  that  the  use  of  more 
advanced  turbulence  model,  for  example  the  two-layer 
turbulence  model,  will  remove  all  the  defeciencies  arised 
trom  using  wall  function  method. 

Comparisions  of  the  predictions  with  the  measured 
data  for  both  axial  and  radial  velocity  components  are 
shown  in  Figs.VII-12  to  VII-15.  The  agreements  between 
measured  data  and  the  predictions  are  fairly  good  except 
the  velocity  components  are  a  little  overpredicted  in  the 
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tail  region  and  in  the  near  wake  region  due  to  the 
neglection  of  pressure  gradient  effect  in  the  turbulence 
model  and  in  the  wall  functions.  However,  the  thickening 
of  boundary  layer  in  the  tail  and  in  the  near  wake  region 
is  well  predicted.  It  is  noted  that  as  shown  in  these 
figures,  the  change  of  forebody  shape  does  not  seriously 
influence  the  boundary  layer  development  in  the  region  of 
middle  body,  thereby,  that  in  the  stern  and  wake  region. 

As  shown  in  Fig.VII-16  to  VII-18,  the  predictions  of 
turbulent  kinetic  energy  agree  fairly  well  with 
experimental  measurements.  However,  when  the 
overprediction  of  velocity  components  in  the  regions  of 
tail  and  near  wake  is  considered,  the  agreements  of 
predictions  with  experimental  data  are  not  improved 
results  over  the  predictions  presented  in  the  chapter  VI. 
If  the  velocity  components  are  accurately  predicted,  it  is 
anticipated  that  the  turbulent  kinetic  energy  in  the 
regions  of  tail  and  near  wake  will  be  overpredicted. 

Overall  agreements  between  measurements  and 
predictions  made  by  the  present  numerical  method  are 
encouraging.  For  the  most  of  body  and  near  wake  region, 
the  pressure  and  velocity  components  are  fairly  well 
predicted.  However,  some  considerations  must  be  made  on 
the  selection  of  turbulence  model  especially  on  the  the 


near  wall  treatment  in  order  to  accuately  predi ct  the 
turbulent  quantities  for  the  flow  involving  pressure 
gradient  and  streamline  curvature  effects .  The  use  of 
more  elaborated  wall  functions  or  introductioii  of  more 
advanced  turbulence  model  is  indeed  needed  to  properly 
simulate  the  flow  field  in  these  situations.  These 
improvements  on  turbulence  model  are  remained  for  future 
study  and  are  not  pursued  further  in  the  present  study. 
However,  the  present  calculation  method  provides  the 
predictions  of  whole  turbulent  flow  field  past 
axisymmetric  bodies  which  are  not  yet  reported  in  the 
literature.  The  present  calculations  also  gives  an 
example  of  the  successful  use  of  two  velocity  staggered 
grid  method  in  the  calculation  of  fluid  flow  in  complex 
geometry  which  has  not  been  reported  in  the  literature 
since  the  original  works  of  Maliska  and  Raithby  [30] . 
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CHAPTER  VIII 

CONCLUSIONS  AND  RECOMMENDATIONS 

A  numerical  study  of  laminar  and  turbulent  flow  past 
finite  two  dimensional  and  axisymmetric  bodies  is 
presented.  The  available  calculation  methods  for 
incompressible  flow  in  a  complex  geometry,  grid  generation 
techniques  and  turbulence  models  are  discussed.  The 
numerical  formulations  and  grid  generation  techniques  and 
solution  procedures  employed  in  the  present  study  are 
presented.  Calculations  are  performed  for  laminar  and 
turbulent  flow  past  a  thin  flat  plate,  turbulent  flow  past 
axisymmetric  bodies  with  changing  body  shapes,  solution 
domains,  turbulence  models  and  numerical  methods.  Major 
contributions  and  findings  of  present  study  may  be 
summarized  as  follows; 

(1)  Wake  function  method  for  the  prediction  of 
turbulent  flow  past  a  flat  plate  is  developed.  An 
accurate  prediction  of  wake  flow  past  a  flat  plate  is  made 
without  the  detailed  calculations  of  the  near  wall  region 
and  the  laminar  wake  region. 


(2)  Numerical  solutions  of  a  complete  flow  field 
past  a  finite  flat  plate  and  axisymmetric  bodies  which 


152 


include  the  leading  edge  interaction,  the  boundary  layer 
development  on  the  body,  the  traling  edge  interaction  and 
the  wake  development,  are  obtained.  The  influence  of 
inlet  boundary  conditions  and  the  location  of  inlet 
solution  boundary  on  the  solution  which  is  observed  in  the 
half  body  calculations  is  removed  in  the  full  body 
calculation.  It  is  observed  that  the  magnitude  of  the 
leading  edge  interaction  is  large,  but  the  interacting 
zone  is  confined  to  a  relatively  small  region.  It  is  also 
found  that  the  partially  parabolic  form  of  Navier-Stokes 
equations  is  not  appropriate  for  the  simulation  of  the 
strong  leading  edge  interaction  although  it  is  suitable 
for  the  simulation  of  the  boundary  layer  on  the  body,  the 
trailing  edge  interaction  and  the  wake  development. 

(3)  Comparisions  of  predictions  by  the  finite 
analytic  method  and  by  the  finite  volume  method  are  made. 
Although  the  finite  analytic  method  gives  better 
convergence  histories  and  more  accurate  predictions,  the 
differences  are  found  to  be  quantitatively  very  small. 

Both  methods  well  predict  the  important  features  of  flow 
field.  It  is  observed  that  the  finite  volume  method 
always  slightly  overpredicts  the  velocity  components  on 
the  boundary  layer  as  well  as  in  the  wake  region. 
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(4)  It  is  found  that  use  of  one  velocity  staggered 
grid  method  with  finite  volume  method  gives  a  limit  in 
convergence.  The  introduction  of  two  velocities  staggered 
method  removes  this  deficiency.  However,  the  two  velocity 
staggered  grid  method  requires  more  storage,  more 
computing  time  and  complicated  programing. 

(5)  The  defect  of  the  k  -  e  turbulence  model  with 
the  Couette  flow  type  wall  function  method  for  the 
prediction  of  transition  of  laminar  flow  to  turbulent 
flow,  turbulent  boundary  layer  involving  pressure  gradient 
and  streamline  curvature  effects  is  addressed.  A  strange 
transition  near  the  leading  edge  was  predicted  and  the 
velocity  components  are  overpredicted  in  the  adverse 
pressure  gradient  boundary  layer  region.  As  compared  with 
the  Couette  flow  type  wall  function  method,  the  two-point 
wall  function  method  well  predicts  the  velocity  components 
in  the  adverse  pressure  gradient  boundary  layer  region  by 
use  of  the  general  law  of  the  wall  in  which  the  pressure 
gradient  effect  on  the  flow  in  the  wall  region  are  taken 
into  consideration.  However,  the  use  of  the  general  law 
of  the  wall  is  limited  and  can  not  be  confidently  applied 
to  the  laminar  flow  region  near  the  leading  edge  and  the 
favorable  pressure  gradient  region.  It  is  observed  that 
the  turbulent  kinetic  energy  are  overpredicted  in  the  tail 


region  of  axisymmetric  bodies .  This  fact  shows  the  the 
defect  k  -  e  turbulence  model  in  the  simulation  of  flow 
field  involving  pressure  gradient  and  streamline  curvature 
effects . 

(6)  Improvements  of  predictions  by  the  introduction 
of  two-layer  turbulence  model  are  observed.  The 
predictions  by  the  two-layer  model  not  only  provide  the 
numerical  solutions  of  viscous-affected  near  wall  region 
without  excessive  computational  efforts,  but  also  the  two- 
layer  model  can  be  applied  to  separated  flow  region  where 
the  wall  function  method  can  not  be  confidently  applied. 
The  use  of  two-layer  model  is  highly  recommended. 

However,  it  is  found  that  the  two-layer  model  slightly 
overpredicts  the  velocity  components  in  the  tail  region  of 
axisymmetric  bodies. 

Although  present  calculation  methods  are  sucessfully 
applied  to  the  simulation  of  various  flow  situations, 
there  exist  much  rooms  for  further  investigations  which 
may  be  summarized  as  follows; 

(1)  The  one  velocity  staggered  grid  method  used  in 
chapter  VI  has  limitations  in  the  general  applications. 

The  two  velocities  staggered  grid  method  employed  in  the 
chapter  VII  removes  the  difficulties  encountered  in  the 
use  of  one  velocity  staggered  method.  However,  this 
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method  requires  a  complicated  programing,  more  storages 
and  a  more  computing  time.  An  introduction  of  newly 
developed  calculation  methods  like  algebraic  manipulation 
method  by  Karki  and  Patankar  [36]  or  the  momentum 
interpolation  method  by  Majumdar  [19]  will  be  alternative 
choices  for  removing  these  deficiencies .  The  introduction 
of  these  methods  are  straight foward.  The  use  of  multi¬ 
grid  method  will  improve  the  convergence. 

(2)  The  origin  of  overprediction  of  velocity 
components  in  the  tail  region  of  axisymmetric  bodies  in 
the  two-layer  calculations  should  be  clearly  investigated 
in  order  to  confidently  apply  this  turbulence  model  to 
general  flow  situations. 

(3)  The  introduction  of  more  advanced  turbulence 
model,  or  at  least  more  elaborated  wall  functions  is 
needed  for  the  proper  prediction  of  complete  turbulent 
flows  past  two  dimensional  or  axisymmetric  bodies.  The 
use  of  two-layer  model  in  this  purpose  may  be  an 
appropriate  choice  although  it  is  questionable  whether 
this  model  may  properly  simulate  the  transition  of  laminar 
flow  to  tur’julent  flow. 

(4)  The  present  calculation  methods  can  be  easily 
applied  to  three  dimensional  situations  with  minor 


modifications . 


Table  1 


Definition  of  Variables 


<1> 

C 

> 

u 

9p  9  9u 

1  9  ,  9v 

r9y^''^ef9y> 

29k 

39x 

V 

^  +^lt 

9p  9  ,  9u, 

1 9  9v 

r9y^^^ef9y^ 

~  ^^cPef^ 

29k 

39y 

k 

G 

-  pe 

e 

+ 

e 

CeiG  k 

e2 

-  P  "=£2  k 

G 


+  2  a. 


oxj  . 


=  0.09,  Ok  =  1.0,  Oe  =  1.3,  Cel  =  1.44,  Ce2  =  1.92. 


a^  =  0,  r  =  1  for  Cartesian  coordinate  system, 
a^  =  1,  y  =  r  for  cylindrical  coordinate  system. 
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Table  2.  Solution  Domains  for  Calculations  in  Chapter-IV 
Type  of  Flow  Re  Grid 

Turbulent  2.48x10®  57x31  -0.6  8.57  1.0 

Re  -  Reynolds  number  based  on  the  free  stream  velocity 
and  body  length . 

Xi  =  location  of  upstream  boundary. 

X(j[  =  location  of  downstream  boundary. 

=  location  of  upper  boundary. 

^The  origin  of  coordinate  system  is  located  at  the  trailing 
edge  of  plate. 


Table  3.  Solution  Domains  for  Calculations  in  Chapter-V® 


Type  of  Flow 

Re 

Grid 

Xd 

Yu 

Laminar 

10® 

135x35 

-1.25 

14.6 

12.7 

Turbulent 

2.48x10® 

120x41 

-1.15 

12.3 

3.9 

^Variables  defined  in  Table  2  and  text. 
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Table  4.  Solution  Domains  for  Calculations  in  Chapter-VI^ 

{  Wall  Function  Method  ) 


Type  of  Body 

Re 

Grid 

^9 

Xd 

Yu 

Afterbody-1 

6. 6x10^ 

96x35 

0.5 

10.9 

1.60 

Afterbody-2 

6.8x10^ 

96x35 

0.6 

10.9 

1.55 

Afterbody-5 

9.3x106 

96x35 

0.6 

10.9 

1.15 

^Variables  defined  in  Table  2  and  text. 

Table  5.  Solution  Domains  for 

(  Two  Layer 

Calculations  in  Chapter-VI® 
Model  ) 

Type  of  Body 

Re 

Grid 

|H 

Xd 

Yu 

Afterbody-1 

6.6x106 

85x53 

0.5 

6.75 

3.64 

Afterbody-2 

6.8x106 

85x53 

0.6 

5.60 

3.71 

Afterbody-3 

6.0x10^ 

90x53 

0.6 

12.4 

4.20 

Afterbody-5 

9.3x10® 

85x53 

0.6 

5.60 

2.73 

^Variables  defined  in  Table  2  and  text. 
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Table  6.  Solution  Domains  for  Calculations  in  Chapter-VII^ 


Type  of  Body 

Re 

Grid 

Xi 

Xd 

Yu 

Afterbody-1 

6. 6x106 

151x36 

-2.23 

9.29 

1 . 77 

Afterbody-2 

6.8x106 

151x36 

-2.23 

9.29 

1.80 

Afterbody -3 

6.0x106 

151x36 

-2.23 

9.29 

2.04 

Afterbody -5 

9.3x10^ 

151x36 

-2.23 

9.29 

1.33 

^Variables  defined  in  Table  2  and  text . 
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y  *  y 

^  -^u 

X  *  X. 

1 

Calculation  Region 

X  -  x^ 

a)  Full  Body  Calculation 


y  -  y. 


X  =•  Xj 


Calculation  Region 


X  “  X, 


b)  Half  Body  Calculation 


Fig.II-1.  Solution  Domains 
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X/L 

a)  Laminar  Flow. 


b)  Turbulent  Flow, 


Fig.II-3.  Numerical  Grids  for  A  Flat  Plate. 
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a)  Afterbody- 1 
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b)  Afterbody-2 

Numerical  Grids  for  Axisymmetric  Bodies 
(  Half  Body  Calculation  ) 
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c)  Afterbody-3 


X/L 


d)  Afterbody-5 


Fig.II-4.  Numerical  Grids  for  Axisymmetric  Bodies 
(  Half  Body  Calculation  ) 
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b)  Afterbody- 2 


Fig.ii-5.  Numerical  Grids  for  Axisymmetric  Bodies 
(  Full  Body  Calculation  ) 
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c)  Afterbody-3 


X/L 


d)  Afterbody-5 


NuTOrlcai  Grids  for  Aaisynmetrl 
(  Full  Body  Calculation  ) 
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Fig. 11-7.  A  Local  Element  for  Finite  Analytic  Method 
in  Boundary-Fitted  Coordinate  System 


Fig.II-8.  A  Local  Element  for  Finite  V< 
in  Boundary-Fitted  Coordinat 
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Fig.IV-l,  Computational  Domain 


Fig,IV-3.  Pressure  Distribution  on  the  Plate 
and  along  the  Wake  Centerline  (  y'*' 


Fig.IVT4.  Pressure  Distribution  Given  in  Ref . [70] 
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Fig.IV-6.  Wake  Centerline  Velocity  Distribution 


Fig.lv-7. 


Distribution 


.  IV-8.  Velocity  Profiles  in  the  Near  wake 


Shear  Stress  Distributions 


127mtn 


Turbulent  Kinetic  Energy  Distributions 
in  the  Near  Wake 


Fig.V“l.  Pressure  Distribution  (  Laminar  Flow  ) 


Fig.V-2.  Pressure  Distribution  near  the  Trailing 
Edge  (  Laminar  Flow  ) 
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Fig.V-4.  Wake  Centerline  Velocity  Distribution 
(  Laminar  Flow  ) 
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Fig.V-7.  Friction  Velocity  Distribution 
(  Turbulent  Flow  ) 


Fig.V-8.  Wake  Centerline  Velocity  Distribution 
(  Turbulent  Flow  ) 
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Fig.VI-l.  Notations  and  Solution  Domains 


Fig. VI-2.  Convergence  History 
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Fig. VI-4.  Pressure  Convergence  History 
(  Finite  Volume  Method  ) 
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Fig. VI-5.  Friction  Velocity  Convergence  History 
(  Finite  Analytic  Method  ) 
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Fig. VI-6.  Friction  Velocity  Convergence  History 
(  Finite  Volume  Method  ) 
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Fig.VI-7.  Wake  Centerline  Velocity  Convergence 
History  (  Finite  Analytic  Method  ) 
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Fig. VI-8.  Wake  Centerline  Velocity  Convergence 
History  (  Finite  Volume  Method  ) 


Fig. VI-10.  Surface  Pressure  Distribution 
(  Afterbody-2  ) 


Fig. VI-11.  Surface  Pressure  Distribution 
(  Afterbody-5  ) 


Fig. VI-12.  Friction  Velocity  Distribution 
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Fig. VI-14.  Friction  Velocity  Distribution 
(  Afterbody-5  ) 
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Fig. VI-24.  Surface  Pressure  Distribution 
(  Afterbody-1  ) 


Fig. VI-25.  Surface  Pressure  Distribution 
(  Afterbody-2  ) 


Fig. VI-26.  Surface  Pressure  Distribution 
(  Afterbody-3  ) 
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Fig. VI-28.  Friction  Velocity  Distribution 
(  Afterbody-1  ) 


Fig. VI-29.  Friction  Velocity  Distribution 
(  Afterbody-2  ) 
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Friction  Velocity  Distribution  for 
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Friction  Velocity  Distribution  for 
Afterbody-5  Given  in  Ref. [56] 


:A=o.9 


X/l=0.7950  X/l=0.9150  X/l=0.93 


XA=0-7036  Vl?0.8308  VL=0.8727 


SDUI(|/(0)|-||) 


KOU1)|/(0)|-<1) 


xoui(|/(OI|-«) 


•a 

•o 

z 

< 

•.> 

•  ' 

L 

1 

>1 

•0 

r® 

» 

0 

;a 

4) 

M 

< 

0) 

V 

•o 

5 

rH 

^  > 

4-1 

“•  * 

0 

Oi 

-• 

>1 

• 

4J 

■r^ 

V 
0 
•— 1 
« 

m 

> 

• 

r* 

9 

n 

•  O 

1 

z  ■ 

H 

> 

o» 

h 

XOUJM/(OH-«) 


x»iij)(/(OV'» 


X/L=0.846  XA=0.934 


O  Ciparlmant 
— Colculotion 
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Fig.VII-7.  Surface  Pressure  Distribution 
(  Afterbody-5  ) 
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Fig.VlI-9.  Friction  Velocity  Distribution 
(  Afterbody-2  ) 


Fig.VII-11.  Friction  Velocity  Distribution 
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APPENDIX 

COMPUTER  CODES 


Wall  function  Method 


1-1)  Geometric  Data  Files 

1-2)  Computer  Program  for  the  Generation  of 
Boundary-Fitted  Coordinates 

1-3)  Computer  Program  for  Flow  Calculaution 
by  Finite  Analytic  Method 

1-4)  Computer  Program  for  Flow  Calculation 
by  Finite  Volume  Method 


2 .  Two-Layer  Model 


2-1)  Geometric  Data  Files 

2-2)  Computer  Program  for  the  Generation  of 
Boundary-Fitted  Coordinates 

2-3)  Computer  Program  for  Flow  Calculation 
by  Finite  Analytic  Method 
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SUBROUTINE  SOURCEU 


SU(I,J)*SS>REFF*(TNUY*(OUy+DVX)  CALL  GEOCOE 

+2.*TMUX«D0X-2./3.*DKX)  DVET-0. 5*(V(I, J+1)-V(I, J-1) ) 

DVXI-0.5*(V(I+1,J)-V(I-1,J) ) 

RETURN  DVX-AJI*(B11*DVXI+B12*DVET) 

end  DVy-AJI*(B21*DVXI4^B22*DVET) 
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Full  Body  Calculation 


1-1)  Computer  Program  for  the  Generation  of 
Body  Shape  and  X-Coordinates 

1-2)  Geometric  Data  Files 

(  Generated  by  Program  1-1  ) 

1-3)  Computer  Program  for  the  Generation  of 
Boundary-Fitted  Coordinates 

1-4)  Computer  Program  for  Input  Data  for 
Flow  calculation 

1-5)  Computer  Program  for  Flow  Calculation 
by  Finite  Volume  Method 
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